diff --git a/.gitignore b/.gitignore index b3e3cc5..25e6771 100644 --- a/.gitignore +++ b/.gitignore @@ -4,4 +4,9 @@ /release/ __pycache__/ orisys.egg-info/ -/.cursor/ \ No newline at end of file +/.cursor/ +/dist_pyd/ +/dist_single_pyd/ +/old_orisys/ +/dist_nuitka/ +/build_nuitka/ \ No newline at end of file diff --git a/Readme.md b/Readme.md index 74812e6..63b40b3 100644 --- a/Readme.md +++ b/Readme.md @@ -1,26 +1,36 @@ -# OrisysSDK API 文档 +# orisys API 文档 -本文件提供了 OrisysSDK 的各类方法 +本文件提供了 `orisys` 套件的中文文档 ## 使用方法 -pip install . +1. 根据平台和python 版本下载对应的`.whl`文件 +2. 运行pip命令安装： +```bash +pip install orisys-{version}-{python-version}-none-{os}_amd64.whl +``` +3. 并测试导入： +```python +import orisys +print(orisys.__version__) +``` --- ## `Sensor` 类别 ### 描述 -创建一个传感器实例 +`orisys`套件的核心类，用于从视频文件或实时摄像头中读取传感器图像，创建一个具体的传感器对象（实例）。并完成形变、光流及力信息的计算与管理。 ### 输入参数 -* **cam_id** (`int | str`, 可选): 传感器 ID、序列号或视频路径。默认为 `0`。 +* **vid_src** (`int | str`, 必填): 传感器 ID、相机编号或视频路径。默认为 `0`。 * **isstitch** (`bool`，可选): 是否启用图像拼接，默认为 `True`。 * **cuda** (`bool`，可选): 是否使用 GPU cuda 加速，默认为 `False`。 +* **verbose** (`bool`，可选): 是否输出初始化及运行过程中的详细日志信息，默认为 `False`。 +* **config_name** (`str`，可选): 内部算法或传感器配置名称，用于加载对应的参数配置，默认为 `prototype`。 ### 返回 -* `Sensor` 对象 +* `Sensor` 实例，用于后续图像采集、形变计算和数据读取操作。 ### 示例 ```python # Example 1： 用相机编号创建 from orisys import Sensor -from src.sensor import Sensor sensor = Sensor(0) # Example 2： 用视频文件创建 @@ -29,63 +39,203 @@ sensor = Sensor("/data/data.mp4") --- -## 2. `get_img` 方法 +## 2. `Sensor.get_img` 方法 ### 描述 -从传感器获取图像，进行复眼图像提取和拼接。 +采集并预处理一帧图像。 + +### 返回 +* `np.ndarray`：处理后的灰度图(`uint8`，尺寸为`config_resolution_stitch`对应分辨率)。 +* `None`：当输入源读取失败（例如视频结束、摄像头断开）时返回。 +### 状态更新 +* 成功时更新`sensor`类里面的灰度图。 +### 示例 + +```python +img = sensor.get_img() +if img is None: + break +``` +--- + +## 3. `Sensor.compute_deformation` 方法 + +### 描述 + +基于当前帧图像计算光流/形变向量场，并通过 HHD 分解将形变场分解为法向/切向分量，同时估计法向力、切向力与扭矩等标量输出。该方法会更新 +`sensor`的形变场与力相关状态，并返回当前帧的光流场。 +调用前置条件：需要先通过`get_img()`成功更新图像。 ### 返回 -* 提取（`issitch = False`）或拼接后（`issitch = True`）的图像 +无 +### 示例 + +```python +img = sensor.get_img() +if img is None: + break + +sensor.compute_deformation() +``` --- +## 4. `Sensor.compute_contact` 方法 -## 3. `compute_deformation` 方法 +### 描述 +基于已计算的法向形变场进行接触区域分析，输出接触质心与轮廓等信息，并生成接触强度图。 + +调用前置条件：需先执行`compute_deformation()`。 + +### 输入参数 +* **min_contour_area** (`int`，可选): 多接触点模式下的最小轮廓面积阈值，默认为 `3000`。 +* **min_magnitude** (`float`，可选): 最小光流幅度阈值，低于该阈值的区域不参与接触分析，默认为 `1.0`。 +* **expansion_threshold** (`float`，可选): 扩张/收缩判定阈值，用于接触区域分割，默认为 `0.01`。 +* **downsample_factor** (`int`，可选): 下采样因子，用于平滑与加速计算，默认为 `1`。 +* **kernel_size** (`tuple`，可选): 形态学操作核大小，默认为 `(5, 5)`。 +* **denoise_sigma** (`float`，可选): 高斯降噪强度；设为` 0 `可禁用降噪，默认为 `15`。 + +### 返回 +`dict`为接触分析结果字典，至少包含以下键： +* **expansion_mask**: 扩张区域掩码（`np.ndarray`）。 +* **contraction_mask**: 收缩区域掩码（`np.ndarray`）。 +* **divergence_map**: 散度图（`np.ndarray`）。 +* **centroid**: 质心坐标（单点模式为 `(x, y)`；多点模式为列表）。 +* **centroid_found**: 是否检测到接触区域（`bool`）。 +* **contour**: 最大轮廓或轮廓列表（可选字段，存在时返回）。 + +### 状态更新 + +* **sensor.centroid**: 接触质心（单点为 `(x, y)`，多点为列表）。 +* **sensor.centroid_found**: 是否找到质心（`bool`）。 + +### 示例 +```python +sensor.compute_contact() +centroid, contour, depth = sensor.read_info( + sensor.info.CENTROID, + sensor.info.CONTOUR, + sensor.info.DEPTH, +) +``` +--- +## 4. `Sensor.read_info` 方法 ### 描述 -计算传感器变形 - ---- -## 4. `read_info` 方法 - -### 描述 - -获取指定类型的传感器数据，返回数量和顺序与输入参数一致。 +按指定的信息类型（`sensor.info`枚举值）批量读取传感器状态数据，并按传入顺序返回对应值的元组。 +该方法通常用于在完成形变计算后，读取光流场、分量场、力、帧率等信息。 ### 输入参数 * **args**: 任意数量的 `Sensor.info` 枚举，用于指定需要获取的数据类型。支持的枚举值及对应数据如下： - * VRAW: Optional[np.ndarray] # 原始光流场, shape=(400, 400, 2) - * VNORMAL: Optional[np.ndarray] # 法向光流场, shape=(400, 400, 2) - * VSHEAR: Optional[np.ndarray] # 切向光流场, shape=(400, 400,2) - * DEPTH: Optional[np.ndarray] # 深度场, shape=(400, 400,1) +### 支持的枚 info keys： +向量/图像类 +* VRAW: 原始光流场, `np.ndarray`，shape=(H, W, 2) +* VNORMAL: 法向光流场, `np.ndarray`，shape=(H, W, 2) +* VSHEAR: 切向光流场, `np.ndarray`，shape=(H, W,2) +* DEPTH: 深度场（可选）, 未启用/未生成时返回空数组`np.array([])`。 - * FNORMAL: Optional[np.float64] # 法相力信息 - * FSHEARX: Optional[np.float64] # 横向切向力分量 - * FSHEARY: Optional[np.float64] # 纵向切向力分量 +标量类 +* FNORMAL: 法向力，`float`。 +* FSHEARX: X切向力，`float`。 +* FSHEARY: Y切向力，`float`。 +* FPS: 输出帧率，`float`。 - * FPS: Optional[float] # 输出帧率 - * TIMESTAMP: Optional[str] # 时间戳 +时间戳 +* TIMESTAMP: 时间戳字符串（HH:MM:SS），`str`。 - * TimeStamp: Optional[float] # 传感器时间戳，单位s +注：H/W 由当前配置分辨率决定。 ### 返回 -* 所请求的传感器数据（返回数量和顺序与输入参数一致） +`Tuple` + +按传入的 `sensor.info` 参数顺序返回一个元组， +返回值数量与有效参数数量一致。 + +- 向量/图像类数据：`np.ndarray` + - 若当前流程中尚未生成该数据，返回空数组 `np.array([])` +- 标量类数据：`float` +- 时间戳：`str`（`TIMESTAMP`，格式为 `HH:MM:SS`） + +### 示例 +```python +# Example 1：读取标量数据 +fps, fn, fx, fy = sensor.read_info( + sensor.info.FPS, + sensor.info.FNORMAL, + sensor.info.FSHEARX, + sensor.info.FSHEARY, +) + +# Example 2： 读取单个向量场（注意解包方式） +flow = sensor.read_info(sensor.info.VRAW) +``` --- -## 5. `disconnect` 方法 +## 5. `Sensor.disconnect` 方法 ### 描述 -断开及释放资源。 +释放传感器相关资源并停止后台异步线程。 调用后将关闭视频/摄像头输入源，并销毁所有 OpenCV 窗口。 + +建议在程序退出前调用以确保资源被正确释放。 --- +--- +# 配置标定相关 +## 1. `list_config` 方法 + +### 描述 +列出内置传感器配置参数字典。 + +### 返回 + +* `list`: 内置传感器配置参数字典的名字列表。 +--- +## 2. `export_config` 方法 + +### 描述 +导出当前传感器实例使用的配置参数字典。 + +### 输入参数 +* **config_name**: 需要导出的配置参数字典名字，`str` +* **path**: 导出的路径，默认`"./config/exported_config.json"` + +### 返回 + +* `dict`: 导出的配置参数字典。 + +### 保存 + +* 导出内置的`config_name`配置参数字典至`path` + +--- --- ## 常见问题解答 (FAQ) -TBC \ No newline at end of file +### Q1:为什么必须先调用`get_img()`，再调用`compute_deformation()`？ +A：`compute_deformation()`依赖最近一次获取并拼接的图像数据进行形变与力计算。 +在调用该方法前，必须先通过`get_img()`更新当前帧图像，否则无法进行有效计算。 + +### Q2:为什么调用了`read_info()`，但读到的数据是空的？ +A：`read_info()`仅用于读取当前已计算的数据状态，不会触发任何计算。 +部分字段依赖特定计算步骤生成，例如： +* `VRAW / VNORMAL / VSHEAR`依赖`compute_deformation()`； +* `DEPTH / CENTROID / CONTOUR`依赖`compute_contact()`。 + +### Q3:`compute_contact()`为一定要调用吗？什么时候需要？ +A：当需要接触相关信息（如接触质心、轮廓或深度/散度图）时，必须调用`compute_contact()`。 +若仅关注形变场或力信息（如光流、法向力、切向力），则无需调用该方法。 + +### Q4:`read_info()`的返回值顺序是如何确定的？ +A：`read_info()`按传入的`sensor.info`参数顺序返回一个元组。 +返回值顺序与参数顺序一一对应，与内部存储顺序无关。 + +### Q5:如果没有标定文件，或拼接结果不正确，会发生什么？ +A：若未指定或未找到标定文件，系统会在需要时创建新的拼接标定数据并保存到默认路径。 +拼接结果异常通常与标定数据、拼接配置或输入图像质量有关，建议重新生成拼接标定并检查配置参数。 \ No newline at end of file diff --git a/build_multi_version.py b/build_multi_version.py index 8a116f6..40a1a9b 100644 --- a/build_multi_version.py +++ b/build_multi_version.py @@ -11,6 +11,9 @@ import platform from pathlib import Path +# Package version +VERSION = "0.2.0" + # Python versions to build for PYTHON_VERSIONS = ['3.8', '3.9', '3.10'] @@ -102,8 +105,10 @@ def obfuscate_code(): # PyArmor 9.x generates version-independent obfuscated code # Output directly to dist_obfuscated/src (PyArmor will create orisys subdirectory) + # Use --obf-code 1 to enable code obfuscation (default is 1, but explicit is better) + # Use --exclude to exclude __init__.py files that don't need obfuscation (optional) if IS_WINDOWS: - cmd = 'pyarmor gen --output dist_obfuscated/src --recursive src/orisys' + cmd = 'pyarmor gen --output dist_obfuscated/src --recursive --obf-code 1 src/orisys' print(f"Running: {cmd}") result = subprocess.run(cmd, capture_output=True, text=True, shell=True) else: @@ -112,6 +117,7 @@ def obfuscate_code(): 'gen', '--output', 'dist_obfuscated/src', '--recursive', + '--obf-code', '1', 'src/orisys' ] print(f"Running: {' '.join(cmd)}") @@ -124,6 +130,26 @@ def obfuscate_code(): print("Code obfuscation successful!") print(result.stdout) + + # Verify obfuscation: check that obfuscated files exist and original .py files are obfuscated + obfuscated_dir = Path('dist_obfuscated/src/orisys') + if obfuscated_dir.exists(): + py_files = list(obfuscated_dir.rglob('*.py')) + print(f"Found {len(py_files)} obfuscated .py files") + # Check a sample file to verify it's obfuscated (should contain PyArmor markers) + if py_files: + sample_file = py_files[0] + try: + with open(sample_file, 'rb') as f: + content = f.read(100) # Read first 100 bytes + # PyArmor obfuscated files typically start with specific markers + if b'__pyarmor__' in content or b'pyarmor' in content.lower(): + print(f"✓ Verified: {sample_file.name} is obfuscated") + else: + print(f"⚠ Warning: {sample_file.name} may not be properly obfuscated") + except Exception as e: + print(f"Could not verify obfuscation: {e}") + return True @@ -132,24 +158,25 @@ def copy_additional_files(): print("\n=== Copying additional files ===") # Ensure configs directory exists in obfuscated output - # PyArmor should copy it, but we verify and copy if needed + # IMPORTANT: Only copy JSON files, NOT .py files (those should be obfuscated by PyArmor) configs_src = Path('src/orisys/configs') configs_dst = Path('dist_obfuscated/src/orisys/configs') if configs_src.exists(): - print(f"Ensuring configs directory is copied...") + print(f"Ensuring configs directory is copied (JSON files only)...") if not configs_dst.exists(): - print(f"Copying configs directory from {configs_src} to {configs_dst}") - shutil.copytree(configs_src, configs_dst, dirs_exist_ok=True) + configs_dst.mkdir(parents=True, exist_ok=True) + + # Only copy JSON files, not any .py files + json_files_src = list(configs_src.glob('*.json')) + json_files_dst = list(configs_dst.glob('*.json')) + + if len(json_files_dst) < len(json_files_src): + print(f"Copying JSON files from configs directory ({len(json_files_src)} files)") + for json_file in json_files_src: + shutil.copy(json_file, configs_dst / json_file.name) else: - # Verify JSON files are present - json_files_src = list(configs_src.glob('*.json')) - json_files_dst = list(configs_dst.glob('*.json')) - if len(json_files_dst) < len(json_files_src): - print(f"Updating configs directory (found {len(json_files_dst)}/{len(json_files_src)} JSON files)") - shutil.copytree(configs_src, configs_dst, dirs_exist_ok=True) - else: - print(f"Configs directory already present with {len(json_files_dst)} JSON files") + print(f"Configs directory already present with {len(json_files_dst)} JSON files") # Copy config directory (if exists at root level) if os.path.exists('config'): @@ -157,11 +184,45 @@ def copy_additional_files(): shutil.copytree('config', 'dist_obfuscated/config', dirs_exist_ok=True) # Copy README and other metadata files - files_to_copy = ['README.md', 'Readme(未完成).md', 'requirements.txt', 'LICENSE'] + files_to_copy = ['README.md', 'requirements.txt', 'LICENSE'] for file in files_to_copy: if os.path.exists(file): print(f"Copying {file}...") shutil.copy(file, f'dist_obfuscated/{file}') + + # Verify that no original .py files are in dist_obfuscated (except obfuscated ones) + print("\n=== Verifying obfuscated output ===") + obfuscated_py_files = list(Path('dist_obfuscated/src/orisys').rglob('*.py')) + print(f"Found {len(obfuscated_py_files)} .py files in obfuscated output") + + # Check if any files look like original source (not obfuscated) + original_like_files = [] + for py_file in obfuscated_py_files: + # Skip __init__.py files as they might be simple + if py_file.name == '__init__.py': + continue + try: + with open(py_file, 'rb') as f: + content = f.read(500) # Read first 500 bytes + # Check for common Python source patterns that shouldn't be in obfuscated code + if b'def ' in content and b'__pyarmor__' not in content: + # Might be original source + if b'import ' in content or b'from ' in content: + # Check if it has readable imports (obfuscated code has different structure) + text_content = content.decode('utf-8', errors='ignore') + if 'import ' in text_content and '__pyarmor__' not in text_content: + original_like_files.append(py_file) + except Exception as e: + print(f"Warning: Could not check {py_file}: {e}") + + if original_like_files: + print(f"⚠️ Warning: Found {len(original_like_files)} files that may not be obfuscated:") + for f in original_like_files[:5]: # Show first 5 + print(f" - {f}") + if len(original_like_files) > 5: + print(f" ... and {len(original_like_files) - 5} more") + else: + print("✓ All .py files appear to be obfuscated") def create_setup_for_obfuscated(python_version=None): @@ -217,7 +278,7 @@ for pkg in packages: setup( name="orisys", - version="0.2.0", + version="{VERSION}", description="Orisys - Optical Tactile Sensor SDK (Protected)", author="Your Name", package_dir={{"": "src"}}, @@ -228,6 +289,8 @@ setup( python_requires="{version_constraint}", # Include PyArmor runtime files zip_safe=False, + # Exclude original source files - only include obfuscated files from dist_obfuscated + exclude_package_data={{'': ['*.pyc', '*.pyo', '__pycache__']}}, options={{ 'bdist_wheel': {{ 'python_tag': '{py_tag}', @@ -243,7 +306,7 @@ setup( def create_manifest_in(): - """Create MANIFEST.in to ensure JSON files in configs are included.""" + """Create MANIFEST.in to ensure JSON files in configs are included and exclude original source.""" print("\n=== Creating MANIFEST.in ===") manifest_content = """# Include all JSON files in configs directory @@ -255,6 +318,12 @@ recursive-include src *.pyd recursive-include src *.so recursive-include src *.dll recursive-include src *.dylib + +# Exclude original source files (only obfuscated .py files should be included) +# The obfuscated .py files from PyArmor will be included automatically +# But we want to make sure no original source files are accidentally included +prune src/orisys +recursive-include src/orisys *.py """ manifest_path = Path('dist_obfuscated/MANIFEST.in') @@ -335,7 +404,7 @@ def build_wheel_for_version(python_version): return False print(f"✅ Wheel for Python {python_version} on {PLATFORM_NAME} built successfully!") - print(f" Expected: orisys-0.2.0-{py_tag}-none-{plat_tag}.whl") + print(f" Expected: orisys-{VERSION}-{py_tag}-none-{plat_tag}.whl") return True @@ -444,7 +513,7 @@ def main(): for wheel in Path(f'dist/{PLATFORM_NAME}').glob('*.whl'): print(f" - {wheel.name}") print("\nInstall with:") - print(f" pip install dist/{PLATFORM_NAME}/orisys-0.2.0-*.whl") + print(f" pip install dist/{PLATFORM_NAME}/orisys-{VERSION}-*.whl") if __name__ == '__main__': diff --git a/build_nuitka.py b/build_nuitka.py new file mode 100644 index 0000000..fca9ee4 --- /dev/null +++ b/build_nuitka.py @@ -0,0 +1,129 @@ +import os +import sys +import subprocess +import shutil +from pathlib import Path + +def install_nuitka(): + print("Checking Nuitka installation...") + try: + import nuitka + print("Nuitka is installed.") + except ImportError: + print("Installing Nuitka...") + subprocess.check_call([sys.executable, "-m", "pip", "install", "nuitka", "ordered-set", "zstandard"]) + +def build(): + install_nuitka() + + version = "0.2.0" + src_dir = Path("src/orisys") + output_dir = Path("build_nuitka") + dist_dir = Path("dist_nuitka") + + # Clean previous builds + if output_dir.exists(): shutil.rmtree(output_dir) + if dist_dir.exists(): shutil.rmtree(dist_dir) + output_dir.mkdir() + dist_dir.mkdir() + + print("\n=== Compiling with Nuitka ===") + # Set PYTHONPATH to include src so nuitka finds the package + env = os.environ.copy() + env["PYTHONPATH"] = "src" + os.pathsep + env.get("PYTHONPATH", "") + + # Nuitka arguments + # --module: compile as module/package + # --include-package=orisys: include the whole package + # --output-dir: where to put result + # src/orisys: the target directory + cmd = [ + sys.executable, "-m", "nuitka", + "--module", + "--include-package=orisys", + "--output-dir=" + str(output_dir), + "--remove-output", + "--no-pyi-file", + "src/orisys" + ] + + print(f"Running: {' '.join(cmd)}") + subprocess.check_call(cmd, env=env) + + # Find the generated .pyd file + # It should be something like orisys.cp310-win_amd64.pyd + pyd_files = list(output_dir.glob("*.pyd")) + if not pyd_files: + raise RuntimeError("Nuitka did not produce a .pyd file!") + + pyd_file = pyd_files[0] + print(f"Generated compiled extension: {pyd_file}") + + print("\n=== Preparing Wheel Structure ===") + # We want the structure: + # orisys/ + # __init__.pyd <-- renamed from orisys.pyd + # configs/ + # *.json + + wheel_root = dist_dir / "wheel_root" + pkg_dir = wheel_root / "orisys" + pkg_dir.mkdir(parents=True) + + # Move and rename pyd + dest_pyd = pkg_dir / "__init__.pyd" + shutil.copy(pyd_file, dest_pyd) + print(f"Copied to {dest_pyd}") + + # Copy JSON configs + configs_src = src_dir / "configs" + configs_dst = pkg_dir / "configs" + configs_dst.mkdir() + + json_count = 0 + for json_file in configs_src.glob("*.json"): + shutil.copy(json_file, configs_dst) + json_count += 1 + print(f"Copied {json_count} JSON config files") + + # Create setup.py for packaging + setup_py_path = dist_dir / "setup.py" + setup_content = f""" +from setuptools import setup, Distribution, find_packages + +class BinaryDistribution(Distribution): + def has_ext_modules(self): + return True + +setup( + name="orisys", + version="{version}", + description="Orisys - Optical Tactile Sensor SDK (Nuitka Compiled)", + author="Orisys", + packages=['orisys'], + package_dir={{'': 'wheel_root'}}, + package_data={{ + 'orisys': ['*.pyd', 'configs/*.json'], + }}, + include_package_data=True, + distclass=BinaryDistribution, + zip_safe=False, +) +""" + with open(setup_py_path, "w") as f: + f.write(setup_content) + + print("\n=== Building Wheel ===") + subprocess.check_call([sys.executable, "setup.py", "bdist_wheel"], cwd=dist_dir) + + # Find the wheel + dist_output = dist_dir / "dist" + wheels = list(dist_output.glob("*.whl")) + if wheels: + print(f"\n[SUCCESS] Wheel created: {wheels[0]}") + else: + print("\n[ERROR] Wheel not found.") + +if __name__ == "__main__": + build() + diff --git a/build_pyd.py b/build_pyd.py new file mode 100644 index 0000000..4a015cd --- /dev/null +++ b/build_pyd.py @@ -0,0 +1,902 @@ +""" +Build script to create wheel distributions with .pyd/.so files instead of .py files. + +This script uses Cython to compile Python source files into extension modules (.pyd on Windows, .so on Linux), +making the source code completely hidden in the final distribution. + +Requirements: + - Cython + - A C compiler (Visual Studio Build Tools on Windows, GCC on Linux) +""" +import os +import shutil +import subprocess +import sys +import platform +from pathlib import Path +import tempfile + + +# Package version +VERSION = "0.2.0" + +# Python versions to build for +PYTHON_VERSIONS = ['3.8', '3.9', '3.10'] + +# Detect platform +IS_WINDOWS = platform.system() == 'Windows' +IS_LINUX = platform.system() == 'Linux' +PLATFORM_NAME = 'win' if IS_WINDOWS else 'linux' if IS_LINUX else 'unknown' + +# Extension for compiled modules +EXT_EXTENSION = '.pyd' if IS_WINDOWS else '.so' + + +def clean_build_dirs(): + """Remove previous build artifacts.""" + dirs_to_clean = ['build', 'dist', 'orisys.egg-info', 'dist_pyd', 'build_pyd'] + for dir_name in dirs_to_clean: + if os.path.exists(dir_name): + print(f"Cleaning {dir_name}...") + try: + # Try to remove with error handling + shutil.rmtree(dir_name, ignore_errors=False, onerror=_handle_remove_error) + except (PermissionError, OSError) as e: + print(f" ⚠️ Warning: Could not fully remove {dir_name}: {e}") + print(f" Some files may be locked. Continuing anyway...") + # Try again with ignore_errors=True + try: + shutil.rmtree(dir_name, ignore_errors=True) + except Exception: + pass + + # Ensure dist directory exists with platform subdirectories + os.makedirs(f'dist/{PLATFORM_NAME}', exist_ok=True) + + +def _handle_remove_error(func, path, exc_info): + """Handle errors during directory removal.""" + # If it's a permission error, try to change permissions and retry + if func is os.rmdir: + try: + # Try to change permissions and remove + os.chmod(path, 0o777) + os.rmdir(path) + except Exception: + # If that fails, just skip it + pass + elif func is os.remove: + try: + os.chmod(path, 0o777) + os.remove(path) + except Exception: + pass + + +def check_cython(): + """Check if Cython is installed.""" + try: + import Cython + print(f"✓ Cython {Cython.__version__} is available") + return True + except ImportError: + print("❌ Cython is not installed!") + print(" Install it with: pip install cython") + return False + + +def check_c_compiler(): + """Check if a C compiler is available.""" + if IS_WINDOWS: + # On Windows, check for cl.exe (Microsoft Visual C++) + # Note: cl.exe might not be in PATH even if installed + # We'll do a basic check but the actual compilation will provide better errors + cl_path = shutil.which('cl.exe') + if cl_path: + print("✓ C compiler (cl.exe) found in PATH") + return True + else: + print("⚠️ C compiler (cl.exe) not found in PATH") + print(" (This is OK if Visual Studio is installed but not in PATH)") + print(" The build will attempt to use the compiler anyway...") + return True # Return True to allow build attempt - it will fail with better error if missing + else: + # On Linux/Mac, check for gcc + try: + result = subprocess.run( + ['gcc', '--version'], capture_output=True, text=True, timeout=5 + ) + if result.returncode == 0: + print(f"✓ C compiler (gcc) is available") + return True + except (FileNotFoundError, subprocess.TimeoutExpired): + pass + + print("❌ C compiler (gcc) not found!") + print("\n Install with:") + if IS_LINUX: + print(" sudo apt-get install build-essential # Ubuntu/Debian") + print(" sudo yum install gcc gcc-c++ # CentOS/RHEL") + else: + print(" Install Xcode Command Line Tools:") + print(" xcode-select --install") + return False + + +def get_python_executable(version): + """Get the Python executable for a specific version based on platform.""" + current_version = f"{sys.version_info.major}.{sys.version_info.minor}" + if current_version == version: + return [sys.executable] + + if IS_WINDOWS: + username = os.environ.get("USERNAME", "") + possible_paths = [ + f'py -{version}', + f'python{version}', + f'C:\\Users\\{username}\\anaconda3\\envs\\py{version.replace(".", "")}\\python.exe', + f'C:\\Users\\{username}\\anaconda3\\envs\\py{version.replace(".", "")}-env\\python.exe', + f'C:\\Users\\{username}\\anaconda3\\envs\\python{version.replace(".", "")}\\python.exe', + f'C:\\Users\\{username}\\anaconda3\\envs\\orisys-py{version.replace(".", "")}\\python.exe', + f'C:\\Users\\{username}\\miniconda3\\envs\\py{version.replace(".", "")}\\python.exe', + f'C:\\Python{version.replace(".", "")}\\python.exe', + f'C:\\Program Files\\Python{version.replace(".", "")}\\python.exe', + ] + else: + possible_paths = [ + f'python{version}', + f'/usr/bin/python{version}', + f'/usr/local/bin/python{version}', + f'~/.pyenv/versions/{version}/bin/python', + f'/opt/python{version}/bin/python', + ] + + for path in possible_paths: + try: + if path.startswith('~'): + path = os.path.expanduser(path) + + if IS_WINDOWS and path.startswith('py -'): + result = subprocess.run( + path.split() + ['-c', 'import sys; print(f"{sys.version_info.major}.{sys.version_info.minor}")'], + capture_output=True, text=True, shell=True + ) + if result.returncode == 0 and result.stdout.strip() == version: + return path.split() + else: + if not os.path.exists(path): + continue + + result = subprocess.run( + [path, '-c', 'import sys; print(f"{sys.version_info.major}.{sys.version_info.minor}")'], + capture_output=True, text=True + ) + if result.returncode == 0 and result.stdout.strip() == version: + return [path] + except (FileNotFoundError, subprocess.SubprocessError, PermissionError): + continue + + return None + + +def create_cython_setup_files(): + """Create setup.py files for Cython compilation.""" + print("\n=== Creating Cython setup files ===") + + # Find all Python files to compile (excluding __init__.py and configs) + src_dir = Path('src/orisys') + py_files = [] + + for py_file in src_dir.rglob('*.py'): + # Skip __init__.py files (they need to stay as .py for package structure) + if py_file.name == '__init__.py': + continue + # Skip getDensity.py (excluded from compilation) + if py_file.name == 'getDensity.py': + continue + # Skip configs directory (JSON files only) + if 'configs' in py_file.parts and py_file.name != '__init__.py': + continue + + # Get relative path from src/orisys + rel_path = py_file.relative_to(src_dir) + py_files.append(rel_path) + + print(f"Found {len(py_files)} Python files to compile") + + # Create dist_pyd directory structure + dist_dir = Path('dist_pyd/src/orisys') + dist_dir.mkdir(parents=True, exist_ok=True) + + # Copy ALL Python source files to dist_pyd (needed for Cython compilation) + # Cython needs the source files to compile them + print("Copying Python source files for compilation...") + for py_file in src_dir.rglob('*.py'): + rel_path = py_file.relative_to(src_dir) + # dst_file should be in dist_pyd/src/orisys/ (not dist_pyd/src/) + dst_file = dist_dir / rel_path + dst_file.parent.mkdir(parents=True, exist_ok=True) + shutil.copy(py_file, dst_file) + if py_file.name == '__init__.py': + print(f" Copied: {rel_path} (will remain as .py)") + else: + print(f" Copied: {rel_path} (will be compiled to .pyd/.so)") + + # Copy configs directory (JSON files) + configs_src = src_dir / 'configs' + configs_dst = dist_dir / 'configs' + if configs_src.exists(): + shutil.copytree(configs_src, configs_dst, dirs_exist_ok=True) + print(f" Copied: configs/ directory") + + # Create setup.py for Cython compilation + # Build the py_files list as a string + py_files_str = '[\n' + for py_file in py_files: + # Convert path to module name (e.g., "sensor.py" -> "orisys.sensor") + module_path = str(py_file).replace('\\', '/').replace('.py', '').replace('/', '.') + file_path = f"src/orisys/{py_file}" + py_files_str += f' ("orisys.{module_path}", "{file_path}"),\n' + py_files_str += ']' + + # Build setup.py content using string concatenation to avoid f-string escaping issues + # Note: We need to handle Windows multiprocessing properly + # The if __name__ == '__main__': guard is required for Windows multiprocessing support + setup_content = '''from setuptools import setup, Extension +from Cython.Build import cythonize +from pathlib import Path +import numpy +import sys +import platform +import multiprocessing +import os + +# Package version +VERSION = "''' + VERSION + '''" + +# Get include directories +numpy_include = numpy.get_include() + +# On Windows, find Python library for explicit linking if needed +if platform.system() == 'Windows': + import os + python_libs_dir = os.path.join(sys.base_prefix, 'libs') + python_version = f"{sys.version_info.major}{sys.version_info.minor}" + python_lib_name = f"python{python_version}" + # Check if library file exists + python_lib_path = os.path.join(python_libs_dir, f"{python_lib_name}.lib") + has_python_lib = os.path.exists(python_lib_path) if python_libs_dir and os.path.exists(python_libs_dir) else False +else: + python_libs_dir = None + python_lib_name = None + python_lib_path = None + has_python_lib = False + +# List of Python files to compile (excluding __init__.py) +py_files = ''' + py_files_str + ''' + +# Create Extension objects +extensions = [] +for module_name, file_path in py_files: + # Setup library directories and libraries for Windows + library_dirs = [] + libraries = [] + extra_link_args = [] + + if platform.system() == 'Windows' and has_python_lib: + # Explicitly link Python library on Windows + library_dirs = [python_libs_dir] + # Use libraries parameter to link Python library + # This is the standard way setuptools expects + libraries = [python_lib_name] + + ext = Extension( + module_name, + [file_path], + include_dirs=[numpy_include], + library_dirs=library_dirs, + libraries=libraries, + # Add any additional compiler flags here + extra_compile_args=[], + extra_link_args=extra_link_args, + ) + extensions.append(ext) + +# Determine number of threads for parallel compilation +# On Windows, disable parallel compilation to avoid multiprocessing spawn issues +# On Linux/Mac, use parallel compilation for faster builds +if platform.system() == 'Windows': + nthreads = 1 # Single-threaded on Windows to avoid multiprocessing issues +else: + nthreads = 4 # Parallel compilation on Linux/Mac + +# Windows multiprocessing requires this guard +if __name__ == '__main__': + # Required for Windows multiprocessing support + multiprocessing.freeze_support() + + # On Windows, ensure 64-bit build if Python is 64-bit + if platform.system() == 'Windows': + import sysconfig + # Check if Python is 64-bit + if sys.maxsize > 2**32: + # Force 64-bit build by setting distutils environment + import distutils.util + # This should force 64-bit build + try: + # Try to set the platform to win-amd64 + import distutils._msvccompiler + # The compiler should automatically detect 64-bit Python + pass + except: + pass + + setup( + name="orisys", + version=VERSION, + description="Orisys - Optical Tactile Sensor SDK (Compiled)", + author="Your Name", + ext_modules=cythonize( + extensions, + compiler_directives={ + 'language_level': "3", + 'embedsignature': True, + # Optimize for size and speed + 'optimize.use_switch': True, + 'optimize.unpack_method_calls': True, + }, + build_dir="build_pyd", + # Create .pyd/.so files instead of .py + nthreads=nthreads, # Parallel compilation (disabled on Windows) + ), + zip_safe=False, + include_package_data=True, + package_data={ + "orisys": ["configs/*.json"], + }, + ) +''' + + setup_file = Path('dist_pyd/setup.py') + with open(setup_file, 'w') as f: + f.write(setup_content) + + print(f"Created setup.py with {len(py_files)} modules to compile") + return py_files + + +def compile_with_cython(python_version): + """Compile Python files to .pyd/.so using Cython.""" + print(f"\n{'='*60}") + print(f"Compiling with Cython for Python {python_version} on {PLATFORM_NAME}") + print('='*60) + + python_cmd = get_python_executable(python_version) + + if not python_cmd: + print(f"⚠️ Python {python_version} not found on {PLATFORM_NAME}, skipping...") + return False + + print(f"Using Python: {' '.join(python_cmd)}") + + # Verify Python version + result = subprocess.run( + python_cmd + ['-c', 'import sys; print(f"Python {sys.version}")'], + capture_output=True, text=True + ) + print(result.stdout.strip()) + + # Check if Cython is installed in this Python + result = subprocess.run( + python_cmd + ['-c', 'import Cython; print(f"Cython {Cython.__version__}")'], + capture_output=True, text=True + ) + if result.returncode != 0: + print(f"❌ Cython not installed for Python {python_version}") + print(f" Install with: {' '.join(python_cmd)} -m pip install cython") + return False + print(result.stdout.strip()) + + # Check if numpy is installed (required for compilation) + result = subprocess.run( + python_cmd + ['-c', 'import numpy; print(f"NumPy {numpy.__version__}")'], + capture_output=True, text=True + ) + if result.returncode != 0: + print(f"❌ NumPy not installed for Python {python_version}") + print(f" Install with: {' '.join(python_cmd)} -m pip install numpy") + return False + print(result.stdout.strip()) + + os.chdir('dist_pyd') + + # On Windows, ensure we use the correct architecture (64-bit) + # Check Python architecture first + if IS_WINDOWS: + arch_check = subprocess.run( + python_cmd + ['-c', 'import sys; print("64bit" if sys.maxsize > 2**32 else "32bit")'], + capture_output=True, text=True + ) + if arch_check.returncode == 0: + python_arch = arch_check.stdout.strip() + print(f"Python architecture: {python_arch}") + if python_arch == "64bit": + # Force 64-bit build on Windows by setting environment variable + # This tells distutils to use the 64-bit toolchain + os.environ['DISTUTILS_USE_SDK'] = '1' + # Also try to find and use vcvars64.bat if available + # But setuptools should handle this automatically + print(" Using 64-bit build configuration") + + # Build extensions + # Don't use --inplace to avoid directory issues + # Build to build/lib directory, then copy to src/orisys + cmd = python_cmd + ['setup.py', 'build_ext'] + print(f"Running: {' '.join(cmd)}") + result = subprocess.run(cmd, capture_output=True, text=True) + + os.chdir('..') + + if result.returncode != 0: + print(f"❌ Error compiling with Cython") + + # Check for common compiler errors + error_output = result.stderr + "\n" + result.stdout + if "cl.exe" in error_output and "failed" in error_output: + print("\n ⚠️ C compiler (cl.exe) error detected!") + print(" This usually means:") + print(" 1. Visual Studio Build Tools are not installed") + print(" 2. The compiler is not in your PATH") + print(" 3. The Visual Studio environment is not set up") + print("\n Solutions:") + print(" - Install Visual Studio Build Tools with C++ support") + print(" - Or open 'Developer Command Prompt for VS' and run from there") + print(" - Or run: vcvars64.bat before building") + + print(f"\n Full error output:\n{result.stderr}") + if result.stdout and result.stderr != result.stdout: + print(f"\n stdout:\n{result.stdout}") + return False + + print(f"✅ Compilation successful for Python {python_version}!") + + # Copy compiled files from build directory to src/orisys + build_lib_dir = Path('build') + target_dir = Path('src/orisys') + + # Find the build directory (it may be build/lib.win-amd64-cpython-310 or similar) + compiled_files = [] + if build_lib_dir.exists(): + # Find the actual build output directory + for lib_dir in build_lib_dir.glob('lib.*'): + orisys_build_dir = lib_dir / 'orisys' + if orisys_build_dir.exists(): + # Copy all .pyd/.so files to src/orisys + print(f" Copying compiled files from {orisys_build_dir} to {target_dir}...") + target_dir.mkdir(parents=True, exist_ok=True) + for compiled_file in orisys_build_dir.glob(f'*{EXT_EXTENSION}'): + # Copy with original filename (e.g., sensor.cp310-win_amd64.pyd) + target_file = target_dir / compiled_file.name + shutil.copy(compiled_file, target_file) + compiled_files.append(target_file) + + # Also create a symlink/copy with base name (e.g., sensor.pyd) + # Python can import from either name, but base name is more reliable + base_name = compiled_file.stem.split('.')[0] # Get 'sensor' from 'sensor.cp310-win_amd64' + base_file = target_dir / f"{base_name}{EXT_EXTENSION}" + if not base_file.exists(): + shutil.copy(compiled_file, base_file) + compiled_files.append(base_file) + print(f" Copied: {compiled_file.name} -> {target_file} and {base_file.name}") + else: + print(f" Copied: {compiled_file.name} -> {target_file}") + break + + # Verify files are in the target location + if not compiled_files: + compiled_files = list(target_dir.rglob(f'*{EXT_EXTENSION}')) + if compiled_files: + print(f" Found {len(compiled_files)} existing {EXT_EXTENSION} files in target directory") + + print(f" Created {len(compiled_files)} {EXT_EXTENSION} files") + + # List some example files + if compiled_files: + print(f" Examples: {', '.join([f.name for f in compiled_files[:3]])}") + if len(compiled_files) > 3: + print(f" ... and {len(compiled_files) - 3} more") + else: + print(" ⚠️ Warning: No compiled files found in expected locations") + print(" This may indicate a problem with the build process") + + # Remove original .py files (except __init__.py) after successful compilation + # This ensures only compiled .pyd/.so files remain + print("\n=== Removing original .py source files (keeping only compiled extensions) ===") + removed_count = 0 + for py_file in Path('dist_pyd/src/orisys').rglob('*.py'): + # Keep __init__.py files + if py_file.name == '__init__.py': + continue + # Check if corresponding .pyd/.so file exists + compiled_file = py_file.with_suffix(EXT_EXTENSION) + if compiled_file.exists(): + try: + py_file.unlink() + removed_count += 1 + print(f" Removed: {py_file.relative_to('dist_pyd')}") + except Exception as e: + print(f" Warning: Could not remove {py_file.name}: {e}") + + if removed_count > 0: + print(f"Removed {removed_count} original .py files (compiled versions exist)") + else: + print("No .py files removed (compiled files may not have been created)") + + # Return True if compilation was successful (returncode == 0) + # Even if file verification had issues, compilation itself succeeded + return True + + +def create_package_setup(python_version=None): + """Create setup.py for packaging the compiled .pyd/.so files.""" + print(f"\n=== Creating package setup.py (Python {python_version or '3.8-3.10'}) ===") + + if python_version: + version_constraint = f'=={python_version}.*' + py_tag = f'py{python_version.replace(".", "")}' + else: + version_constraint = '>=3.8,<3.11' + py_tag = 'py38.py39.py310' + + # Before packaging, ensure .pyd files are in src/orisys + # They might still be in build directory, so copy them if needed + print("=== Ensuring .pyd files are in src/orisys before packaging ===") + orisys_dir = Path('dist_pyd/src/orisys') + build_dir = Path('dist_pyd/build') + + # Check if files are in build directory and copy them + copied_count = 0 + if build_dir.exists(): + for lib_dir in build_dir.glob('lib.*'): + orisys_build_dir = lib_dir / 'orisys' + if orisys_build_dir.exists(): + print(f" Found .pyd files in {orisys_build_dir}, copying to src/orisys...") + orisys_dir.mkdir(parents=True, exist_ok=True) + for compiled_file in orisys_build_dir.glob(f'*{EXT_EXTENSION}'): + target_file = orisys_dir / compiled_file.name + # Always copy to ensure we have the latest version + shutil.copy(compiled_file, target_file) + copied_count += 1 + print(f" Copied: {compiled_file.name} -> {target_file}") + + if copied_count > 0: + print(f" Successfully copied {copied_count} .pyd files to src/orisys") + else: + print(" ⚠️ No .pyd files found in build directory to copy") + + # Before packaging, ensure all .py files (except __init__.py) are removed + print("=== Final cleanup: Removing any remaining .py source files ===") + removed_count = 0 + if orisys_dir.exists(): + for py_file in orisys_dir.rglob('*.py'): + # Keep __init__.py files + if py_file.name == '__init__.py': + continue + # Check if corresponding .pyd/.so file exists + compiled_file = py_file.with_suffix(EXT_EXTENSION) + # Check for any .pyd file with this base name + base_name = py_file.stem + has_compiled = False + if compiled_file.exists(): + has_compiled = True + else: + # Check for versioned .pyd files (e.g., sensor.cp310-win_amd64.pyd) + for pyd_file in orisys_dir.glob(f'{base_name}.*{EXT_EXTENSION}'): + has_compiled = True + break + + if has_compiled: + try: + py_file.unlink() + removed_count += 1 + print(f" Removed: {py_file.relative_to('dist_pyd')}") + except Exception as e: + print(f" Warning: Could not remove {py_file.name}: {e}") + else: + # If no compiled version exists, this is a problem + print(f" ⚠️ Warning: {py_file.name} has no compiled version, removing anyway") + try: + py_file.unlink() + removed_count += 1 + except Exception as e: + print(f" Warning: Could not remove {py_file.name}: {e}") + + if removed_count > 0: + print(f"Removed {removed_count} .py files before packaging") + + # Find all compiled .pyd/.so files in src/orisys + compiled_files = list(orisys_dir.rglob(f'*{EXT_EXTENSION}')) if orisys_dir.exists() else [] + print(f"Found {len(compiled_files)} compiled {EXT_EXTENSION} files in src/orisys") + + if len(compiled_files) == 0: + print(" ⚠️ WARNING: No .pyd files found in src/orisys!") + print(" This will cause the wheel to be missing compiled modules.") + + # Create MANIFEST.in to ensure .pyd files are included + manifest_content = f"""# Include all compiled extensions +include src/orisys/*{EXT_EXTENSION} +recursive-include src/orisys *{EXT_EXTENSION} + +# Include JSON config files +recursive-include src/orisys/configs *.json + +# Exclude all .py files except __init__.py +global-exclude *.py +include src/orisys/__init__.py +recursive-include src/orisys/configs __init__.py +""" + manifest_file = Path('dist_pyd/MANIFEST.in') + with open(manifest_file, 'w') as f: + f.write(manifest_content) + print("Created MANIFEST.in to ensure .pyd files are included") + + setup_content = f'''from setuptools import setup, find_packages +from setuptools.command.build_py import build_py as _build_py +from pathlib import Path +import os +import shutil +import glob + +# Read requirements +with open('requirements.txt', 'r') as f: + requirements = [line.strip() for line in f if line.strip() and not line.startswith('#')] + +# Find all packages +packages = find_packages(where="src") + +# Collect package data (JSON files and .pyd/.so files) +# IMPORTANT: .pyd/.so files with version tags (e.g., sensor.cp310-win_amd64.pyd) +# are NOT automatically included by setuptools, so we must explicitly list them +package_data = {{}} +for pkg in packages: + pkg_path = Path('src') / pkg.replace('.', '/') + if pkg_path.exists(): + data_files = [] + + # Explicitly include compiled extensions (.pyd/.so) using glob pattern + # This ensures all .pyd files are included regardless of version tags + for pattern in ['*{EXT_EXTENSION}']: + for ext_file in pkg_path.glob(pattern): + rel_path = ext_file.relative_to(pkg_path) + data_files.append(str(rel_path)) + # Also check subdirectories + for ext_file in pkg_path.rglob(pattern): + rel_path = ext_file.relative_to(pkg_path) + if str(rel_path) not in data_files: + data_files.append(str(rel_path)) + + # Include JSON files from configs directory + if pkg == 'orisys.configs': + for json_file in pkg_path.glob('*.json'): + rel_path = json_file.relative_to(pkg_path) + data_files.append(str(rel_path)) + + # For orisys package, include configs subdirectory JSON files + if pkg == 'orisys': + configs_dir = pkg_path / 'configs' + if configs_dir.exists() and configs_dir.is_dir(): + for json_file in configs_dir.glob('*.json'): + rel_path = json_file.relative_to(pkg_path) + data_files.append(str(rel_path)) + + if data_files: + package_data[pkg] = data_files + print(f"Package {{pkg}} package_data: {{len(data_files)}} files") + +# Custom build_py to ensure .pyd files are copied +class build_py(_build_py): + def run(self): + """Override to ensure .pyd files are included""" + super().run() + # After building, verify .pyd files are in the build directory + for pkg in packages: + pkg_path = Path('src') / pkg.replace('.', '/') + build_pkg_path = Path(self.build_lib) / pkg.replace('.', '/') + if pkg_path.exists() and build_pkg_path.exists(): + # Copy any .pyd files that might have been missed + for ext_file in pkg_path.glob('*{EXT_EXTENSION}'): + build_ext_file = build_pkg_path / ext_file.name + if not build_ext_file.exists(): + shutil.copy(ext_file, build_ext_file) + print(f" Copied {{ext_file.name}} to build directory") + +setup( + name="orisys", + version="{VERSION}", + description="Orisys - Optical Tactile Sensor SDK (Compiled)", + author="Your Name", + package_dir={{"": "src"}}, + packages=packages, + package_data=package_data, + include_package_data=True, + install_requires=requirements, + python_requires="{version_constraint}", + zip_safe=False, + cmdclass={{"build_py": build_py}}, + options={{ + 'bdist_wheel': {{ + 'python_tag': '{py_tag}', + }} + }}, +) +''' + + with open('dist_pyd/setup.py', 'w') as f: + f.write(setup_content) + + print("Package setup.py created successfully!") + + +def build_wheel_for_version(python_version): + """Build the wheel distribution for a specific Python version.""" + print(f"\n{'='*60}") + print(f"Building wheel for Python {python_version} on {PLATFORM_NAME}") + print('='*60) + + python_cmd = get_python_executable(python_version) + + if not python_cmd: + print(f"⚠️ Python {python_version} not found, skipping...") + return False + + # Get platform tag + plat_tag = get_platform_tag() + py_tag = f'py{python_version.replace(".", "")}' + + os.chdir('dist_pyd') + + # Build the wheel + cmd = python_cmd + ['setup.py', 'bdist_wheel', + '--python-tag', py_tag, + '--plat-name', plat_tag, + '--dist-dir', f'../dist/{PLATFORM_NAME}'] + print(f"Running: {' '.join(cmd)}") + result = subprocess.run(cmd, capture_output=True, text=True) + + os.chdir('..') + + if result.returncode != 0: + print(f"❌ Error building wheel:\n{result.stderr}") + if result.stdout: + print(f"stdout:\n{result.stdout}") + return False + + print(f"✅ Wheel for Python {python_version} built successfully!") + print(f" Expected: orisys-{VERSION}-{py_tag}-none-{plat_tag}.whl") + return True + + +def get_platform_tag(): + """Get the platform tag for wheel naming.""" + system = platform.system() + machine = platform.machine().lower() + + if system == 'Windows': + if machine in ('amd64', 'x86_64'): + return 'win_amd64' + elif machine in ('x86', 'i386', 'i686'): + return 'win32' + elif machine in ('arm64', 'aarch64'): + return 'win_arm64' + elif system == 'Linux': + if machine in ('x86_64', 'amd64'): + return 'manylinux2014_x86_64' + elif machine in ('aarch64', 'arm64'): + return 'manylinux2014_aarch64' + elif machine in ('i386', 'i686'): + return 'manylinux2014_i686' + elif system == 'Darwin': + if machine == 'arm64': + return 'macosx_11_0_arm64' + else: + return 'macosx_10_9_x86_64' + + return 'any' + + +def copy_additional_files(): + """Copy non-Python files.""" + print("\n=== Copying additional files ===") + + # Copy requirements.txt + if os.path.exists('requirements.txt'): + shutil.copy('requirements.txt', 'dist_pyd/requirements.txt') + print("Copied requirements.txt") + + # Copy README and LICENSE if they exist + for file in ['README.md', 'LICENSE']: + if os.path.exists(file): + shutil.copy(file, f'dist_pyd/{file}') + print(f"Copied {file}") + + +def main(): + """Main build process.""" + import argparse + + parser = argparse.ArgumentParser(description='Build orisys package with Cython (.pyd/.so files)') + parser.add_argument('--versions', nargs='+', default=PYTHON_VERSIONS, + help='Python versions to build for') + + args = parser.parse_args() + + current_py_version = f"{sys.version_info.major}.{sys.version_info.minor}" + + print("=" * 60) + print("Building Orisys Package with Cython (.pyd/.so files)") + print(f"Platform: {PLATFORM_NAME} ({platform.system()})") + print(f"Current Python: {current_py_version} ({sys.executable})") + print(f"Python Versions: {', '.join(args.versions)}") + print("=" * 60) + + # Check Cython + if not check_cython(): + print("\n❌ Build failed: Cython is required!") + print(" Install with: pip install cython") + sys.exit(1) + + # Check C compiler + if not check_c_compiler(): + print("\n❌ Build failed: C compiler is required!") + print(" Please install the required build tools (see above)") + sys.exit(1) + + # Ensure we're in the right directory + script_dir = Path(__file__).parent + os.chdir(script_dir) + + # Step 1: Clean previous builds + clean_build_dirs() + + # Step 2: Create Cython setup files + py_files = create_cython_setup_files() + if not py_files: + print("\n❌ No Python files found to compile!") + sys.exit(1) + + # Step 3: Copy additional files + copy_additional_files() + + # Step 4: Compile for each Python version + success_count = 0 + for version in args.versions: + # Compile with Cython + if compile_with_cython(version): + # Create package setup.py + create_package_setup(version) + # Build wheel + if build_wheel_for_version(version): + success_count += 1 + + if success_count == 0: + print("\n❌ Build failed - no wheels created!") + sys.exit(1) + + print(f"\n✅ Successfully built {success_count}/{len(args.versions)} wheels") + print("\n" + "=" * 60) + print("✅ Build completed successfully!") + print("=" * 60) + print(f"\nPlatform: {PLATFORM_NAME}") + print(f"Compiled source: dist_pyd/") + print(f"Distribution packages: dist/{PLATFORM_NAME}/") + print("\nGenerated wheels:") + for wheel in Path(f'dist/{PLATFORM_NAME}').glob('*.whl'): + print(f" - {wheel.name}") + print("\nInstall with:") + print(f" pip install dist/{PLATFORM_NAME}/orisys-{VERSION}-*.whl") + print("\nNote: The wheel contains .pyd/.so files instead of .py files!") + + +if __name__ == '__main__': + main() + diff --git a/example1_basic_pipeline.py b/example1_basic_pipeline.py index b14ce84..b6aeb68 100644 --- a/example1_basic_pipeline.py +++ b/example1_basic_pipeline.py @@ -19,11 +19,10 @@ def main(): # 步骤1：创建传感器对象 (必要步骤)，config_name为配置文件名（内置配置文件名）或路径（自定义配置文件路径） # 方式1: 从视频文件读取 - sensor = orisys.Sensor("./stereo2.mp4", isstitch=True, cuda=True, verbose=True,config_name="prototype") + sensor = orisys.Sensor("./s1.mp4", isstitch=True, cuda=True, verbose=True,config_name="prototype") # 方式2: 从摄像头读取（取消注释使用） # sensor = orisys.Sensor(0, cuda=True, verbose=True, config_name="prototype") - sensor.export_config(path="./config/exported_config.json") # 导出配置文件供查看和修改 print("\n按 'q' 退出程序，按 'r' 重置光流追踪器") diff --git a/src/orisys/HHD.py b/src/orisys/HHD.py new file mode 100644 index 0000000..50a3e7e --- /dev/null +++ b/src/orisys/HHD.py @@ -0,0 +1,1356 @@ +''' +Copyright (c) 2015, Harsh Bhatia (bhatia4@llnl.gov) +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: +1. Redistributions of source code must retain the above copyright notice, this +list of conditions and the following disclaimer. +2. Redistributions in binary form must reproduce the above copyright notice, +this list of conditions and the following disclaimer in the documentation +and/or other materials provided with the distribution. +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE +FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +''' + +import numpy +# from .nPoisson import PoissonSolver +# from .StructuredGrid import StructuredGrid +import timeit +import sys +from scipy import signal, spatial +# from .UnstructuredGrid import UnstructuredGrid +# import matplotlib.pyplot as plt + +# Try to import CuPy for GPU acceleration +try: + import cupy as cp + CUPY_AVAILABLE = True + + # Try to import CuPy streams for parallel execution + try: + from cupy.cuda import Stream + STREAMS_AVAILABLE = True + except ImportError: + Stream = None + STREAMS_AVAILABLE = False + + # Create fused elementwise kernel for torque computation + # Computes cross products and combines them in a single kernel pass + try: + _torque_kernel = cp.ElementwiseKernel( + 'float32 r_x, float32 r_y, float32 grid_0, float32 grid_1, ' + 'float32 max_0, float32 max_1, float32 min_0, float32 min_1', + 'float32 torque', + ''' + float cross_p = (grid_0 - max_0) * r_y - (grid_1 - max_1) * r_x; + float cross_n = (grid_0 - min_0) * r_y - (grid_1 - min_1) * r_x; + torque = cross_p + cross_n; + ''', + 'torque_kernel' + ) + KERNEL_AVAILABLE = True + except Exception: + _torque_kernel = None + KERNEL_AVAILABLE = False + +except ImportError: + cp = None + CUPY_AVAILABLE = False + STREAMS_AVAILABLE = False + KERNEL_AVAILABLE = False + Stream = None + _torque_kernel = None + +# --------------------------------------------------------------------------- +class nHHD(object): + + # constructor + def __init__(self, **kwargs): + + args = kwargs.keys() + + if ('grid' in args) == ('points' in args): + raise SyntaxError("nHHD Solver needs either shape of a regular grid or the points in an rectilinear/unstructured grid") + + if ('grid' in args) != ('spacings' in args): + raise SyntaxError("nHHD needs spacings for regular grid") + + if 'grid' in args: + + self.mesh = StructuredGrid(grid=kwargs['grid'], spacings=kwargs['spacings']) + self.psolver = PoissonSolver(grid=self.mesh.dims, spacings=self.mesh.dx) + + # elif 'points' in args: + + # if 'simplices' in args: + # self.mesh = UnstructuredGrid(vertices=kwargs['points'], simplices=kwargs['simplices']) + # else: + # self.mesh = UnstructuredGrid(vertices=kwargs['points']) + + # self.psolver = PoissonSolver(points=self.mesh.vertices, pvolumes=self.mesh.pvolumes) + + self.dim = self.mesh.dim + if (self.dim != 2) and (self.dim != 3): + raise ValueError("nHHD Solver works only for 2 and 3 dimensions") + + self.psolver.prepare(False) + + # create the 3-component decomposition + def decompose(self, vfield, verbose=0): + + if vfield.shape[-1] != self.dim: + raise ValueError("nHHD.decompose requires a valid-dimensional vector field") + + if verbose > 0: + print( '\n ==> nhhd.decompose', vfield.shape) + + if self.dim == 2: + + # compute div curl + (self.div, self.curlw) = self.mesh.divcurl(vfield) + + if verbose > 0: + print( 'div =', self.div.shape, self.div.min(), self.div.max()) + print( 'curlw =', self.curlw.shape, self.curlw.min(), self.curlw.max()) + + # compute potentials + self.nD = self.psolver.solve_gpu(self.div, verbose > 0) #1) + if verbose > 0: + print( 'nD =', self.nD.shape, self.nD.min(), self.nD.max()) + + self.nRu = self.psolver.solve_gpu(self.curlw, verbose > 0) #1) + if verbose > 0: + print( 'nR =', self.nRu.shape, self.nRu.min(), self.nRu.max()) + + # compute fields as gradients of potentials + self.d = self.mesh.gradient(self.nD) + if verbose > 0: + print( 'd =', self.d.shape, self.d.min(), self.d.max()) + + self.r = self.mesh.rotated_gradient(self.nRu) + if verbose > 0: + print( 'r =', self.r.shape, self.r.min(), self.r.max()) + else: + + # compute div curl + (self.div, self.curlu, self.curlv, self.curlw) = self.mesh.divcurl(vfield) + + if verbose > 0: + print( 'div =', self.div.shape, self.div.min(), self.div.max()) + print( 'curlu =', self.curlu.shape, self.curlu.min(), self.curlu.max()) + print( 'curlv =', self.curlv.shape, self.curlv.min(), self.curlv.max()) + print( 'curlw =', self.curlw.shape, self.curlw.min(), self.curlw.max()) + + # compute potentials + self.nD = self.psolver.solve_gpu(self.div, verbose > 0) #1) + if verbose > 0: + print( 'nD =', self.nD.shape, self.nD.min(), self.nD.max()) + + self.nRu = self.psolver.solve_gpu(self.curlu, verbose > 0) + self.nRu *= -1.0 + if verbose > 0: + print( 'nRu =', self.nRu.shape, self.nRu.min(), self.nRu.max()) + + self.nRv = self.psolver.solve_gpu(self.curlv, verbose > 0) + self.nRv *= -1.0 + if verbose > 0: + print( 'nRv =', self.nRv.shape, self.nRv.min(), self.nRv.max()) + + self.nRw = self.psolver.solve_gpu(self.curlw, verbose > 0) + self.nRw *= -1.0 + + if verbose > 0: + print( 'nRw =', self.nRw.shape, self.nRw.min(), self.nRw.max()) + + # compute divergent field as gradient of D potential + self.d = self.mesh.gradient(self.nD) + if verbose > 0: + print( 'd =', self.d.shape, self.d.min(), self.d.max()) + + # temporary use of r + self.r = numpy.stack((self.nRu, self.nRv, self.nRw), axis = -1) + + # compute rotational field as curl of R potential + (cu, cv, cw) = self.mesh.curl3D(self.r) + self.r = numpy.stack((cu, cv, cw), axis = -1) + if verbose > 0: + print( 'r =', self.r.shape, self.r.min(), self.r.max()) + + self.h = numpy.add(self.d, self.r) + numpy.subtract(vfield, self.h, self.h) + + if verbose > 0: + print( 'h =', self.h.shape, self.h.min(), self.h.max()) + + def decompose_gpu(self, vfield, verbose=0): + """ + GPU-optimized decomposition that keeps intermediate results on GPU. + Returns GPU arrays to avoid CPU-GPU transfers. + + Args: + vfield: Vector field (numpy array or CuPy array) + verbose: Verbosity level (default: 0) + + Returns: + Tuple of GPU arrays: (d_gpu, h_gpu, r_gpu, nRu_gpu) for 2D fields + """ + if not CUPY_AVAILABLE: + if verbose > 0: + print("CuPy not available, falling back to CPU decompose") + self.decompose(vfield, verbose) + return None + + if vfield.shape[-1] != self.dim: + raise ValueError("nHHD.decompose requires a valid-dimensional vector field") + + if verbose > 0: + print('\n ==> nhhd.decompose_gpu', vfield.shape) + + # Convert input to GPU if needed + if isinstance(vfield, numpy.ndarray): + vfield_gpu = cp.asarray(vfield, dtype=cp.float32) + else: + vfield_gpu = vfield + + if self.dim == 2: + # Compute div curl on GPU + dx0, dx1 = self.mesh.dx[0], self.mesh.dx[1] + + dudy, dudx = cp.gradient(vfield_gpu[:,:,0], dx0, dx1) + dvdy, dvdx = cp.gradient(vfield_gpu[:,:,1], dx0, dx1) + + div_gpu = dudx + dvdy + curlw_gpu = dvdx - dudy + + if verbose > 0: + print('div =', div_gpu.shape, float(div_gpu.min()), float(div_gpu.max())) + print('curlw =', curlw_gpu.shape, float(curlw_gpu.min()), float(curlw_gpu.max())) + + # Compute potentials on GPU (solve_gpu returns numpy, so we convert back) + nD_cpu = self.psolver.solve_gpu(cp.asnumpy(div_gpu), verbose > 0) + nRu_cpu = self.psolver.solve_gpu(cp.asnumpy(curlw_gpu), verbose > 0) + + # Keep potentials on GPU + nD_gpu = cp.asarray(nD_cpu, dtype=cp.float32) + nRu_gpu = cp.asarray(nRu_cpu, dtype=cp.float32) + + if verbose > 0: + print('nD =', nD_gpu.shape, float(nD_gpu.min()), float(nD_gpu.max())) + print('nR =', nRu_gpu.shape, float(nRu_gpu.min()), float(nRu_gpu.max())) + + # Compute gradients on GPU + ddy, ddx = cp.gradient(nD_gpu, dx0, dx1) + d_gpu = cp.stack((ddx, ddy), axis=-1) + + if verbose > 0: + print('d =', d_gpu.shape, float(d_gpu.min()), float(d_gpu.max())) + + # Compute rotated gradient on GPU + rdy, rdx = cp.gradient(nRu_gpu, dx0, dx1) + rdy *= -1.0 + r_gpu = cp.stack((rdy, rdx), axis=-1) + + if verbose > 0: + print('r =', r_gpu.shape, float(r_gpu.min()), float(r_gpu.max())) + + # Compute h on GPU + h_gpu = d_gpu + r_gpu + h_gpu = vfield_gpu - h_gpu + + if verbose > 0: + print('h =', h_gpu.shape, float(h_gpu.min()), float(h_gpu.max())) + + # Store GPU arrays + self._d_gpu = d_gpu + self._h_gpu = h_gpu + self._r_gpu = r_gpu + self._nRu_gpu = nRu_gpu + + # Also store CPU copies for compatibility + self.d = cp.asnumpy(d_gpu) + self.h = cp.asnumpy(h_gpu) + self.r = cp.asnumpy(r_gpu) + self.nRu = cp.asnumpy(nRu_gpu) + self.nD = nD_cpu + + return (d_gpu, h_gpu, r_gpu, nRu_gpu) + + else: + # For 3D, fall back to CPU for now + if verbose > 0: + print("3D GPU decompose not yet implemented, using CPU") + self.decompose(vfield, verbose) + return None + + def get_force(self, vfield, verbose=False): + + mgrid = numpy.mgrid[0: vfield.shape[0], 0: vfield.shape[1]] + grid = numpy.stack((mgrid[0], mgrid[1]), axis =2) + # print(grid.shape) + + + self.decompose(vfield, verbose) + + self.Sn = numpy.sum(numpy.absolute(self.d)) + # print("Sn is:", self.Sn) + + # self.Stx = abs(numpy.sum(self.h[:, :, 0])) + # self.Sty = abs(numpy.sum(self.h[:, :, 1])) + self.Stx = numpy.sum(self.h[:, :, 0]) + self.Sty = numpy.sum(self.h[:, :, 1]) + # print("Stx is:", self.Stx) + # print("Sty is:", self.Sty) + + nR_max_loc = numpy.unravel_index(numpy.argmax(self.nRu, axis = None), self.nRu.shape) + nR_min_loc = numpy.unravel_index(numpy.argmin(self.nRu, axis = None), self.nRu.shape) + + # print(nR_max_loc) + + r_p = grid - nR_max_loc + r_n = grid - nR_min_loc + + self.Stor = numpy.sum(numpy.cross(r_p, self.r)) + numpy.sum(numpy.cross(r_n, self.r)) + # print("Stor is:", self.Stor) + + output = numpy.array([self.Sn, self.Stx, self.Sty, self.Stor]) + + flag = False + if self.Sn >=300 and self.Stx >= 300 and self.Sty >= 300 and numpy.abs(self.Stor) >= 1000: + flag = True + # print(flag) + return flag,output + # return self.d, self.h, numpy.cross(r_p, self.r) , numpy.cross(r_n, self.r) + + def _init_force_gpu_cache(self, shape): + """ + Initialize cached GPU arrays for force computation. + Call once when shape is known to avoid repeated allocations. + + Args: + shape: Tuple (height, width) of the vector field + """ + if not CUPY_AVAILABLE: + return + + # Cache grid arrays (only need to create once per shape) + self._force_gpu_cache = { + 'shape': shape, + 'mgrid_0': cp.arange(shape[0], dtype=cp.float32)[:, None], + 'mgrid_1': cp.arange(shape[1], dtype=cp.float32)[None, :], + } + # Pre-broadcast for efficiency + self._force_gpu_cache['mgrid_0_bc'] = cp.broadcast_to( + self._force_gpu_cache['mgrid_0'], shape + ).copy() # .copy() to make contiguous + self._force_gpu_cache['mgrid_1_bc'] = cp.broadcast_to( + self._force_gpu_cache['mgrid_1'], shape + ).copy() + + def get_force_gpu_v2(self, vfield, verbose=False): + """ + Further optimized GPU force computation using CuPy ElementwiseKernel. + Fuses all operations into a single kernel for maximum performance. + + Args: + vfield: Vector field (numpy array) + verbose: Print debug information (default: False) + + Returns: + tuple: (flag, output) where: + - flag: Boolean indicating if force thresholds are met + - output: Array [Sn, Stx, Sty, Stor] containing force components + """ + if not CUPY_AVAILABLE: + if verbose: + print("CuPy not available, falling back to CPU get_force") + return self.get_force(vfield, verbose) + + # Decompose the vector field + self.decompose(vfield, verbose) + + # Keep data on GPU - avoid transfers by computing directly + d_gpu = cp.asarray(self.d, dtype=cp.float32) + h_gpu = cp.asarray(self.h, dtype=cp.float32) + r_gpu = cp.asarray(self.r, dtype=cp.float32) + nRu_gpu = cp.asarray(self.nRu, dtype=cp.float32) + + # Use reduction kernels for better performance + # Compute Sn using L1 norm (sum of absolute values) + Sn = float(cp.linalg.norm(d_gpu.ravel(), ord=1)) + + # Compute Stx, Sty + Stx = float(cp.sum(h_gpu[:, :, 0])) + Sty = float(cp.sum(h_gpu[:, :, 1])) + + # Find extrema indices + nR_max_idx = int(cp.argmax(nRu_gpu)) + nR_min_idx = int(cp.argmin(nRu_gpu)) + + shape = nRu_gpu.shape + max_0, max_1 = nR_max_idx // shape[1], nR_max_idx % shape[1] + min_0, min_1 = nR_min_idx // shape[1], nR_min_idx % shape[1] + + # Initialize cache if needed + if not hasattr(self, '_force_gpu_cache') or self._force_gpu_cache is None or self._force_gpu_cache['shape'] != shape: + self._init_force_gpu_cache(shape) + + mgrid_0 = self._force_gpu_cache['mgrid_0_bc'] + mgrid_1 = self._force_gpu_cache['mgrid_1_bc'] + + r_x = r_gpu[:, :, 0] + r_y = r_gpu[:, :, 1] + + # Compute torque with fused sum + Stor = float(cp.sum( + (mgrid_0 - max_0) * r_y - (mgrid_1 - max_1) * r_x + + (mgrid_0 - min_0) * r_y - (mgrid_1 - min_1) * r_x + )) + + self.Sn, self.Stx, self.Sty, self.Stor = Sn, Stx, Sty, Stor + output = numpy.array([Sn, Stx, Sty, Stor], dtype=numpy.float32) + + flag = (Sn >= 300 and Stx >= 300 and Sty >= 300 and abs(Stor) >= 1000) + + return flag, output + + def clear_force_gpu_cache(self): + """Clear the GPU cache for force computation to free memory.""" + if hasattr(self, '_force_gpu_cache'): + self._force_gpu_cache = None + if CUPY_AVAILABLE: + cp.get_default_memory_pool().free_all_blocks() + + def get_shear(self): + + return self.h + + def get_normal(self): + + return self.d + +# --------------------------------------------------------------------------- + +class Timer(object): + + def display(self): + + tseconds = self.etime - self.stime + + print(' [[ elapsed time:', end = "") + if tseconds < numpy.power(10.,-6.): + print("%.3f micro-sec." % (tseconds*numpy.power(10,6)),end = "") + + elif tseconds < numpy.power(10.,-3.): + print("%.3f milli-sec." % (tseconds*numpy.power(10.,3.)),end = "") + + elif tseconds < 60.0: + print("%.3f sec." % (tseconds),end = "") + + else: + m = int(tseconds/ 60.0) + print("%d min. %.3f sec." % (m, (tseconds - 60*m)),end = "") + + print(']]') + + def start(self): + self.stime = timeit.default_timer() + + def end(self, print_time=True): + self.etime = timeit.default_timer() + + if print_time: + self.display() + + def __init__(self): + self.start() + +# --------------------------------------------------------------------------- + +# Try to import CuPy for GPU acceleration +try: + import cupy as cp + CUPY_AVAILABLE = True +except ImportError: + cp = None + CUPY_AVAILABLE = False + +# ------------------------------------------------------------------------------ +class StructuredGrid(object): + + # -------------------------------------------------------------------------- + # -------------------------------------------------------------------------- + def __init__(self, **kwargs): + + ''' + kwargs: + grid: ndarray of grid dimensions (Y,X) or (Z,Y,X) + spacings: ndarray of grid spacings (dy, dx) or (dz, dy, dx) + verbose: verbosity level + ''' + + args = kwargs.keys() + + if ('grid' not in args) or ('spacings' not in args): + raise SyntaxError("Dimensions and spacings of the grid are required") + + verbose = 0 + if 'verbose' in args: + verbose = kwargs['verbose'] + + self.dims = kwargs['grid'] + self.dx = kwargs['spacings'] + self.dim = len(self.dims) + + if self.dim != 2 and self.dim != 3: + raise ValueError("StructuredGrid works for 2D and 3D only") + + if self.dim != len(self.dx): + raise ValueError("Dimensions of spacings should match that of the grid") + + # Initialize GPU cache + self._gpu_divcurl_cache = None + + if verbose > 0: + print(' Initializing', self.dim, 'D structured grid...', end ='') + sys.stdout.flush() + mtimer = Timer() + print(' Done!', end ='') + mtimer.end() + + # -------------------------------------------------------------------------- + # -------------------------------------------------------------------------- + def divcurl(self, vfield, verbose=False): + # compute div and curl of vfield (vector field) + #if (vfield.shape[-1] != self.dim) or (vfield.shape[0:self.dim] - self.dims).any(): + if (vfield.shape[-1] != self.dim) or (vfield.shape[0:self.dim] != self.dims): + raise ValueError("Dimensions of vector field should match that of the grid") + + if verbose: + print(' Computing div and curl...', end ='') + sys.stdout.flush() + mtimer = Timer() + + if self.dim == 2: + + # self.dx = (dy,dx) + dudy, dudx = numpy.gradient(vfield[:,:,0], self.dx[0], self.dx[1]) + dvdy, dvdx = numpy.gradient(vfield[:,:,1], self.dx[0], self.dx[1]) + + numpy.add(dudx, dvdy, dudx) + numpy.subtract(dvdx, dudy, dvdx) + + if verbose: + print(' Done!', end ='') + mtimer.end() + + return (dudx, dvdx) + + elif self.dim == 3: + + # self.dx = (dz,dy,dx) + dudz, dudy, dudx = numpy.gradient(vfield[:,:,:,0], self.dx[0], self.dx[1], self.dx[2]) + dvdz, dvdy, dvdx = numpy.gradient(vfield[:,:,:,1], self.dx[0], self.dx[1], self.dx[2]) + dwdz, dwdy, dwdx = numpy.gradient(vfield[:,:,:,2], self.dx[0], self.dx[1], self.dx[2]) + + numpy.add(dudx, dvdy, dudx) + numpy.add(dudx, dvdz, dudx) + + numpy.subtract(dwdy, dvdz, dwdy) + numpy.subtract(dudz, dwdx, dudz) + numpy.subtract(dvdx, dudy, dvdx) + + if verbose: + print(' Done!', end ='') + mtimer.end() + + return (dudx, dwdy, dudz, dvdx) + + def divcurl_gpu(self, vfield, verbose=False): + """ + GPU-accelerated computation of divergence and curl of vector field. + Uses CuPy for faster gradient computation on GPU. + """ + if not CUPY_AVAILABLE: + if verbose: + print('CuPy not available, falling back to CPU') + return self.divcurl(vfield, verbose) + + if (vfield.shape[-1] != self.dim) or (vfield.shape[0:self.dim] != self.dims): + raise ValueError("Dimensions of vector field should match that of the grid") + + if verbose: + print(' Computing div and curl on GPU...', end ='') + sys.stdout.flush() + mtimer = Timer() + + if self.dim == 2: + # Transfer to GPU + vfield_gpu = cp.asarray(vfield, dtype=cp.float32) + + # Compute gradients on GPU + dudy, dudx = cp.gradient(vfield_gpu[:,:,0], self.dx[0], self.dx[1]) + dvdy, dvdx = cp.gradient(vfield_gpu[:,:,1], self.dx[0], self.dx[1]) + + # Compute div and curl on GPU + div_gpu = dudx + dvdy + curl_gpu = dvdx - dudy + + # Transfer back to CPU + div = cp.asnumpy(div_gpu) + curl = cp.asnumpy(curl_gpu) + + if verbose: + print(' Done!', end ='') + mtimer.end() + + return (div, curl) + + elif self.dim == 3: + # Transfer to GPU + vfield_gpu = cp.asarray(vfield, dtype=cp.float32) + + # Compute gradients on GPU + dudz, dudy, dudx = cp.gradient(vfield_gpu[:,:,:,0], self.dx[0], self.dx[1], self.dx[2]) + dvdz, dvdy, dvdx = cp.gradient(vfield_gpu[:,:,:,1], self.dx[0], self.dx[1], self.dx[2]) + dwdz, dwdy, dwdx = cp.gradient(vfield_gpu[:,:,:,2], self.dx[0], self.dx[1], self.dx[2]) + + # Compute div and curl components on GPU + div_gpu = dudx + dvdy + dwdz + curlu_gpu = dwdy - dvdz + curlv_gpu = dudz - dwdx + curlw_gpu = dvdx - dudy + + # Transfer back to CPU + div = cp.asnumpy(div_gpu) + curlu = cp.asnumpy(curlu_gpu) + curlv = cp.asnumpy(curlv_gpu) + curlw = cp.asnumpy(curlw_gpu) + + if verbose: + print(' Done!', end ='') + mtimer.end() + + return (div, curlu, curlv, curlw) + + def curl3D(self, vfield, verbose=False): + + #if (vfield.shape[-1] != self.dim) or (vfield.shape[0:self.dim] - self.dims).any(): + if (vfield.shape[-1] != self.dim) or (vfield.shape[0:self.dim] != self.dims): + raise ValueError("Dimensions of vector field should match that of the grid") + + if self.dim != 3: + raise ValueError("curl3D works only for 2D") + + if verbose: + print(' Computing curl...', end ='') + sys.stdout.flush() + mtimer = Timer() + + # self.dx = (dz,dy,dx) + dudz, dudy, dudx = numpy.gradient(vfield[:,:,:,0], self.dx[0], self.dx[1], self.dx[2]) + dvdz, dvdy, dvdx = numpy.gradient(vfield[:,:,:,1], self.dx[0], self.dx[1], self.dx[2]) + dwdz, dwdy, dwdx = numpy.gradient(vfield[:,:,:,2], self.dx[0], self.dx[1], self.dx[2]) + + numpy.subtract(dwdy, dvdz, dwdy) + numpy.subtract(dudz, dwdx, dudz) + numpy.subtract(dvdx, dudy, dvdx) + + if verbose: + print(' Done!', end ='') + mtimer.end() + + return (dwdy, dudz, dvdx) + + def curl3D_gpu(self, vfield, verbose=False): + """ + GPU-accelerated computation of curl for 3D vector field. + Uses CuPy for faster gradient computation on GPU. + """ + if not CUPY_AVAILABLE: + if verbose: + print('CuPy not available, falling back to CPU') + return self.curl3D(vfield, verbose) + + if (vfield.shape[-1] != self.dim) or (vfield.shape[0:self.dim] != self.dims): + raise ValueError("Dimensions of vector field should match that of the grid") + + if self.dim != 3: + raise ValueError("curl3D works only for 3D") + + if verbose: + print(' Computing curl on GPU...', end ='') + sys.stdout.flush() + mtimer = Timer() + + # Transfer to GPU + vfield_gpu = cp.asarray(vfield, dtype=cp.float32) + + # Compute gradients on GPU + dudz, dudy, dudx = cp.gradient(vfield_gpu[:,:,:,0], self.dx[0], self.dx[1], self.dx[2]) + dvdz, dvdy, dvdx = cp.gradient(vfield_gpu[:,:,:,1], self.dx[0], self.dx[1], self.dx[2]) + dwdz, dwdy, dwdx = cp.gradient(vfield_gpu[:,:,:,2], self.dx[0], self.dx[1], self.dx[2]) + + # Compute curl components on GPU + curl_u = dwdy - dvdz + curl_v = dudz - dwdx + curl_w = dvdx - dudy + + # Transfer back to CPU + curlu = cp.asnumpy(curl_u) + curlv = cp.asnumpy(curl_v) + curlw = cp.asnumpy(curl_w) + + if verbose: + print(' Done!', end ='') + mtimer.end() + + return (curlu, curlv, curlw) + + def rotated_gradient(self, sfield, verbose=False): + + if (sfield.shape != self.dims): + #if (sfield.shape - self.dims).any(): + raise ValueError("Dimensions of scalar field should match that of the grid") + + if self.dim != 2: + raise ValueError("rotated_gradient works only for 2D") + + if verbose: + print(' Computing rotated gradient...', end ='') + sys.stdout.flush() + mtimer = Timer() + + ddy, ddx = numpy.gradient(sfield, self.dx[0], self.dx[1]) + ddy *= -1.0 + + grad = numpy.stack((ddy, ddx), axis=-1) + + if verbose: + print(' Done!', end ='') + mtimer.end() + + return grad + + def rotated_gradient_gpu(self, sfield, verbose=False): + """ + GPU-accelerated computation of rotated gradient for 2D scalar field. + Uses CuPy for faster gradient computation on GPU. + """ + if not CUPY_AVAILABLE: + if verbose: + print('CuPy not available, falling back to CPU') + return self.rotated_gradient(sfield, verbose) + + if (sfield.shape != self.dims): + raise ValueError("Dimensions of scalar field should match that of the grid") + + if self.dim != 2: + raise ValueError("rotated_gradient works only for 2D") + + if verbose: + print(' Computing rotated gradient on GPU...', end ='') + sys.stdout.flush() + mtimer = Timer() + + # Transfer to GPU + sfield_gpu = cp.asarray(sfield, dtype=cp.float32) + + # Compute gradients on GPU + ddy, ddx = cp.gradient(sfield_gpu, self.dx[0], self.dx[1]) + ddy *= -1.0 + + # Stack and transfer back to CPU + grad_gpu = cp.stack((ddy, ddx), axis=-1) + grad = cp.asnumpy(grad_gpu) + + if verbose: + print(' Done!', end ='') + mtimer.end() + + return grad + + def gradient(self, sfield, verbose=False): + + if (sfield.shape != self.dims): + #if (sfield.shape - self.dims).any(): + raise ValueError("Dimensions of scalar field should match that of the grid") + + if verbose: + print(' Computing gradient...', end ='') + sys.stdout.flush() + mtimer = Timer() + + if self.dim == 2: + + # self.dx = (dy,dx) + ddy, ddx = numpy.gradient(sfield, self.dx[0], self.dx[1]) + grad = numpy.stack((ddx, ddy), axis = -1) + + elif self.dim == 3: + + # self.dx = (dz,dy,dx) + ddz, ddy, ddx = numpy.gradient(sfield, self.dx[0], self.dx[1], self.dx[2]) + grad = numpy.stack((ddx, ddy, ddz), axis = -1) + + if verbose: + print(' Done!', end ='') + mtimer.end() + + return grad + + def gradient_gpu(self, sfield, verbose=False): + """ + GPU-accelerated computation of gradient for scalar field. + Uses CuPy for faster gradient computation on GPU. + """ + if not CUPY_AVAILABLE: + if verbose: + print('CuPy not available, falling back to CPU') + return self.gradient(sfield, verbose) + + if (sfield.shape != self.dims): + raise ValueError("Dimensions of scalar field should match that of the grid") + + if verbose: + print(' Computing gradient on GPU...', end ='') + sys.stdout.flush() + mtimer = Timer() + + # Transfer to GPU + sfield_gpu = cp.asarray(sfield, dtype=cp.float32) + + if self.dim == 2: + # Compute gradients on GPU + ddy, ddx = cp.gradient(sfield_gpu, self.dx[0], self.dx[1]) + grad_gpu = cp.stack((ddx, ddy), axis=-1) + + elif self.dim == 3: + # Compute gradients on GPU + ddz, ddy, ddx = cp.gradient(sfield_gpu, self.dx[0], self.dx[1], self.dx[2]) + grad_gpu = cp.stack((ddx, ddy, ddz), axis=-1) + + # Transfer back to CPU + grad = cp.asnumpy(grad_gpu) + + if verbose: + print(' Done!', end ='') + mtimer.end() + + return grad + + def clear_gpu_cache(self): + """Clear GPU cache to free GPU memory""" + if self._gpu_divcurl_cache is not None: + self._gpu_divcurl_cache = None + if CUPY_AVAILABLE: + cp.get_default_memory_pool().free_all_blocks() + + # -------------------------------------------------------------------------- + # -------------------------------------------------------------------------- + + +import numpy +from scipy import signal, spatial +# from .timer import Timer + +# Try to import CuPy for GPU acceleration +try: + import cupy as cp + import cupyx.scipy.signal as cp_signal + CUPY_AVAILABLE = True +except ImportError: + cp = None + cp_signal = None + CUPY_AVAILABLE = False + + +''' +# ------------------------------------------------------------------------------ + Poisson Solver for 1D, 2D, and 3D scalar fields. + Data may be defined on a rectlinear grid or unstructured points. + + 1) create an object + PoissonSolver(args) + 'solver' : type of solver + 'S' for spatial convolution + only option for unstructured points + 'F' for convolution in frequency space + default option for regular grids + + --> for regular grid + 'grid' : shape of the rectlinear grid + tuple of dimensions (Z,Y,X) + 'spacings' : spacings of the rectlinear grid + tuple of spacings (dz,dy,dx) + + --> for unstrucutred points + 'points' : unstructured points + numpy array of shape (N, dim) + 'pvolumes' : point volumes for each point + numpy array of shape (N, 1) + + 2) call prepare() + this initializes internal datastructures depending upon the mesh. + + Once initialized, solve() can be called any number of times to + solve the Poisson equation on the same domain. If the domain + changes (e.g., a Lagrangian mesh), create a new object. + + 3) call solve(f), where, f is a numpy array defining a scalar field + + for rectlinear grid, f.shape = grid.shape + for unstructured points, f.shape = (points.shape[0], 1) + +# ------------------------------------------------------------------------------ +''' + +class PoissonSolver(object): + + # constructor + def __init__(self, **kwargs): + + args = kwargs.keys() + + if ('grid' in args) == ('points' in args): + raise SyntaxError("Poisson Solver needs either shape of a regular grid or the points in an rectilinear/unstructured grid") + + if ('grid' in args) != ('spacings' in args): + raise SyntaxError("Poisson Solver needs spacings for regular grid") + + if ('points' in args) != ('pvolumes' in args): + raise SyntaxError("Poisson Solver needs point volumes for unstructured points") + + # (default) data type + self.dtype = numpy.float32 + if 'dtype' in args: + self.dtype = kwargs['dtype'] + + if 'grid' in args: + self.ptype = 'G' + self.stype = 'F' # default + self.gdims = kwargs['grid'] + self.gdx = kwargs['spacings'] + self.dim = len(self.gdims) + if self.dim != len(self.gdx): + raise ValueError("Dimensions of spacings should match that of the grid") + + elif 'points' in args: + self.ptype = 'P' + self.stype = 'S' # only option + self.points = kwargs['points'] + self.pvolumes = kwargs['pvolumes'] + self.dim = len(self.points.shape) + if self.pvolumes.shape[0] != self.points.shape[0]: + raise ValueError("Number of pvolumes should match that of the points") + + if 'solver' in args: + if (kwargs['solver'] != 'S' and kwargs['solver'] != 'F'): + raise SyntaxError('Solver can be only of type S(patial) or F(requency)') + + if (kwargs['solver'] == 'F') and ('points' in args): + raise SyntaxError('Only Spatial solver supported for unstructured points') + + if (kwargs['solver'] == 'S') and ('grid' in args): + raise SyntaxError('Frequncy solver strongly suggested for structured grids') + + self.stype = kwargs['solver'] + + if (self.dim > 3): + raise ValueError("Poisson Solver works only for 1, 2, or 3 dimensions") + + self.ready = False + self._gpu_cache = None # Cache for GPU arrays + # if (self.ptype == 'G'): + # print( 'PoissonSolver:', self.dim, 'D grid =', self.gdims, 'with spacings', self.gdx) + # elif (self.ptype == 'P'): + # print( 'PoissonSolver', self.dim, 'D points =', self.points.shape) + + # -------------------------------------------------------------------------- + # create the Green's function in appropriate representation + # -------------------------------------------------------------------------- + + # generate a grid -(N-1) to (N-1) to define distance function + # where N is the size of the original grid + def _create_radialGrid(self, verbose): + + mtimer = Timer() + if verbose: + print( ' - creating distance kernel:', end ='') + + rdims = tuple(2*d - 1 for d in self.gdims) + + R = numpy.indices(rdims, dtype=self.dtype) + for d in range(self.dim): + R[d] += 1-self.gdims[d] + R[d] *= self.gdx[d] + + if self.dim == 1: + numpy.absolute(R[0], R[0]) + R = R[0] + + elif self.dim == 2: + numpy.hypot(R[1], R[0], R[0]) + R = R[0] + + # half of the average grid spacing + zval = 0.5 * (self.gdx[0]+self.gdx[1]) / 2.0 + R[self.gdims[0]-1, self.gdims[1]-1] = zval + + elif self.dim == 3: + numpy.hypot(R[2], R[1], R[1]) + numpy.hypot(R[1], R[0], R[0]) + R = R[0] + + # half of the average grid spacing + zval = 0.5 * (self.gdx[0]+self.gdx[1]+self.gdx[2]) / 3.0 + R[self.gdims[0]-1, self.gdims[1]-1, self.gdims[2]-1] = zval + + if verbose: + print( R.shape, R.min(), R.max()) + mtimer.end() + + return R + + # -------------------------------------------------------------------------- + # scale a function by point volumes for correct integration + def _scale(self, func): + + if self.ptype == 'P': + return numpy.multiply(func, self.pvolumes) + + elif self.ptype == 'G': + sfunc = func * numpy.prod(self.gdx) + + ''' + # volume of the boundary elements is only half the voxel + if self.dim == 1: + sfunc[ 0] *= 0.5 + sfunc[-1] *= 0.5 + elif self.dim == 2: + sfunc[ 0, :] *= 0.5 + sfunc[-1, :] *= 0.5 + sfunc[ :, 0] *= 0.5 + sfunc[ :,-1] *= 0.5 + elif self.dim == 3: + sfunc[ 0, :, :] *= 0.5 + sfunc[-1, :, :] *= 0.5 + sfunc[ :, 0, :] *= 0.5 + sfunc[ :,-1, :] *= 0.5 + sfunc[ :, :, 0] *= 0.5 + sfunc[ :, :,-1] *= 0.5 + ''' + return sfunc + + def prepare(self, verbose=False): + + # ---------------------------------------------------------------------- + gtimer = Timer() + if verbose: + print( '\nInitializing Poisson solver, type =', self.stype) + + # ---------------------------------------------------------------------- + # generate a grid to store pairwise distances to store Green's function + self.G = None # definition depends upon solver type + + # ---------------------------------------------------------------------- + # for Fourier solution, + if self.stype == 'F': + + if self.ptype == 'G': + self.G = self._create_radialGrid(verbose) + + elif self.ptype == 'P': + raise SyntaxError('Frequncy solver not supported for unstructured points') + #self.G = create_distArray() + + # ---------------------------------------------------------------------- + # for spatial solution, + # use the pairwise distance matrix, N x N + # where N is the number of points + elif self.stype == 'S': + + if self.ptype == 'G': + raise SyntaxError('Spatial solver not suggested for structured grids') + #self.G = self.create_points(verbose) + + elif self.ptype == 'P': + # for 1D points, add an extra axis + if (self.dim == 1) and (len(self.points.shape) == 1): + self.points = self.points[:,numpy.newaxis] + + # create pairwise distance matrix now + ltimer = Timer() + self.G = spatial.distance.cdist(self.points, self.points) + + # TODO: use a good approximation here + zval = numpy.power(10.0, -10.0) + self.G += zval * numpy.identity(self.G.shape[0]) + + if verbose: + print(' - created pairwise distance matrix:', self.G.shape, self.G.min(), self.G.max(), end = "") + ltimer.end() + + # ---------------------------------------------------------------------- + # compute the Green's function + ltimer = Timer() + if verbose: + print( ' - computing the Green\'s function:') + + if self.dim == 1: + numpy.multiply(self.G, 0.5, self.G) + + elif self.dim == 2: + numpy.log(self.G, self.G) + numpy.multiply(self.G, (0.5 / numpy.pi), self.G) + + elif self.dim == 3: + numpy.reciprocal(self.G, self.G) + numpy.multiply(self.G, (-0.25 / numpy.pi), self.G) + + if verbose: + print( self.G.shape, self.G.min(), self.G.max(), end = "") + ltimer.end() + + # ---------------------------------------------------------------------- + self.ready = True + + if verbose: + print( 'Poisson solver initialized', end ='') + gtimer.end() + + # -------------------------------------------------------------------------- + # compute the integral solution + # -------------------------------------------------------------------------- + def solve(self, f, verbose=False): + + if not self.ready: + self.prepare(verbose) + + fshape = f.shape + + gtimer = Timer() + if verbose: + print( '\nSolving Poisson Eq.') + + # ---------------------------------------------------------------------- + # regular grid + if self.ptype == 'G': + + #if (self.gdims - fshape).any(): + if self.gdims != fshape: + print( self.gdims, fshape) + raise ValueError("Shape of function should match shape of grid") + + # convolution in frequency domain + if self.stype == 'F': + p = signal.fftconvolve(f, self.G, mode='same') + numpy.multiply(p, numpy.prod(self.gdx), p) + + # convolution in spatial domain + elif self.stype == 'S': + raise SyntaxError('Spatial solver not suggested for structured grids') + #p = (self.G).dot(f.reshape(-1)).reshape(fshape) + + # ---------------------------------------------------------------------- + # unstructured points + elif self.ptype == 'P': + + if (self.points.shape[0] != fshape[0]) or (1 != len(fshape)): + raise ValueError("Shape of function should be (N,1) with N as the number of points") + + # convolution in spatial domain + # Dec 13, 2017 -- pre scaling works for unstructured grids! + if self.stype == 'S': + p = numpy.multiply(f, self.pvolumes) + numpy.dot(self.G, p, p) + #numpy.multiply(p, self.pvolumes, p) + + elif self.stype == 'F': + raise SyntaxError('Frequncy solver not supported for unstructured points') + + # ---------------------------------------------------------------------- + if verbose: + print( 'Poisson solver finished') + gtimer.end() + + return p + + # -------------------------------------------------------------------------- + # GPU-accelerated solve using CuPy + # -------------------------------------------------------------------------- + def solve_gpu(self, f, verbose=False): + """ + GPU-accelerated Poisson solver using CuPy. + Falls back to CPU if CuPy is not available. + + Args: + f: Input field (numpy array) + verbose: Print timing information + + Returns: + p: Solution (numpy array on CPU) + """ + if not CUPY_AVAILABLE: + if verbose: + print("CuPy not available, falling back to CPU solver") + return self.solve(f, verbose) + + if not self.ready: + self.prepare(verbose) + + fshape = f.shape + + gtimer = Timer() + if verbose: + print('\nSolving Poisson Eq. on GPU') + + # ---------------------------------------------------------------------- + # regular grid + if self.ptype == 'G': + if self.gdims != fshape: + print(self.gdims, fshape) + raise ValueError("Shape of function should match shape of grid") + + # convolution in frequency domain using GPU + if self.stype == 'F': + # Initialize GPU cache if needed + if self._gpu_cache is None: + if verbose: + print(' - Transferring Green\'s function to GPU...') + self._gpu_cache = { + 'G_gpu': cp.asarray(self.G, dtype=cp.float32), + 'gdx_prod': numpy.prod(self.gdx) + } + + # Transfer input to GPU + f_gpu = cp.asarray(f, dtype=cp.float32) + + # Perform FFT convolution on GPU + # CuPy's fftconvolve is GPU-accelerated + p_gpu = cp_signal.fftconvolve(f_gpu, self._gpu_cache['G_gpu'], mode='same') + + # Scale by grid spacing (on GPU) + p_gpu *= self._gpu_cache['gdx_prod'] + + # Transfer result back to CPU + p = cp.asnumpy(p_gpu) + + elif self.stype == 'S': + raise SyntaxError('Spatial solver not suggested for structured grids') + + # ---------------------------------------------------------------------- + # unstructured points - not optimized for GPU + elif self.ptype == 'P': + if verbose: + print(' - GPU solver not optimized for unstructured points, using CPU') + return self.solve(f, verbose) + + # ---------------------------------------------------------------------- + if verbose: + print('Poisson solver (GPU) finished') + gtimer.end() + + return p + + def solve_gpu_batch(self, f_list, verbose=False): + """ + Batch GPU-accelerated Poisson solver for multiple inputs. + Solves multiple Poisson equations in a single GPU call. + + Args: + f_list: List of input fields (numpy arrays) + verbose: Print timing information + + Returns: + List of solutions (numpy arrays on CPU) + """ + if not CUPY_AVAILABLE: + if verbose: + print("CuPy not available, falling back to CPU solver") + return [self.solve(f, verbose) for f in f_list] + + if not self.ready: + self.prepare(verbose) + + gtimer = Timer() + if verbose: + print(f'\nSolving {len(f_list)} Poisson Eq. on GPU (batched)') + + # Verify all inputs have same shape + fshape = f_list[0].shape + for f in f_list: + if f.shape != fshape: + raise ValueError("All input fields must have the same shape for batch processing") + + # ---------------------------------------------------------------------- + # regular grid + if self.ptype == 'G': + if self.gdims != fshape: + raise ValueError("Shape of function should match shape of grid") + + # convolution in frequency domain using GPU + if self.stype == 'F': + # Initialize GPU cache if needed + if self._gpu_cache is None: + if verbose: + print(' - Transferring Green\'s function to GPU...') + self._gpu_cache = { + 'G_gpu': cp.asarray(self.G, dtype=cp.float32), + 'gdx_prod': numpy.prod(self.gdx) + } + + results = [] + + # Process all inputs in batch + for i, f in enumerate(f_list): + # Transfer input to GPU + f_gpu = cp.asarray(f, dtype=cp.float32) + + # Perform FFT convolution on GPU + p_gpu = cp_signal.fftconvolve(f_gpu, self._gpu_cache['G_gpu'], mode='same') + + # Scale by grid spacing (on GPU) + p_gpu *= self._gpu_cache['gdx_prod'] + + # Transfer result back to CPU + results.append(cp.asnumpy(p_gpu)) + + if verbose: + print(f' - Solved {i+1}/{len(f_list)}') + + if verbose: + print('Poisson solver (GPU batch) finished') + gtimer.end() + + return results + + elif self.stype == 'S': + raise SyntaxError('Spatial solver not suggested for structured grids') + + # ---------------------------------------------------------------------- + # unstructured points - not optimized for GPU + elif self.ptype == 'P': + if verbose: + print(' - GPU solver not optimized for unstructured points, using CPU') + return [self.solve(f, verbose) for f in f_list] + + # If we reach here, it's an unsupported case + raise ValueError(f"Unsupported point type '{self.ptype}' for GPU batch solver") + + def clear_gpu_cache(self): + """Clear GPU cache to free GPU memory""" + if self._gpu_cache is not None: + self._gpu_cache = None + if CUPY_AVAILABLE: + cp.get_default_memory_pool().free_all_blocks() + + \ No newline at end of file diff --git a/src/orisys/MultiSensor.py b/src/orisys/MultiSensor.py deleted file mode 100644 index 4f0ae2d..0000000 --- a/src/orisys/MultiSensor.py +++ /dev/null @@ -1,315 +0,0 @@ -""" -MultiSensor - Multi-process sensor management for parallel processing. - -This module provides the MultiSensor class for running multiple Sensor -instances in parallel processes, enabling true parallel processing with -independent CUDA contexts. -""" - -import os - - -def _sensor_worker(sensor_id, config, data_queue, stop_event): - """ - Worker function that runs in a separate PROCESS for each sensor. - Each process has its own Python interpreter and CUDA context. - """ - os.environ["OPENCV_VIDEOIO_MSMF_ENABLE_HW_TRANSFORMS"] = "0" - - from .sensor import Sensor - - print(f"[Sensor {sensor_id}] Process started (PID: {os.getpid()})") - - sensor = Sensor( - vid_src=config.get('vid_src'), - isstitch=config.get('isstitch', True), - cuda=config.get('cuda', True), - verbose=config.get('verbose', False), - config_name=config.get('config_name', 'prototype'), - cal_path=config.get('cal_path', None) - ) - compute_contact = config.get('compute_contact', False) - - frame_count = 0 - - def queue_put(data): - try: - while not data_queue.empty(): - try: data_queue.get_nowait() - except: break - data_queue.put_nowait(data) - except: pass - - while not stop_event.is_set(): - try: - # Sync mode - sensor.get_img() - sensor.compute_deformation() - if compute_contact: - sensor.compute_contact() - - # Use read_info to get all sensor data - - img, flow, normal, shear, fps, fnormal, fshearx, fsheary = sensor.read_info( - sensor.info.IMG, - sensor.info.VRAW, - sensor.info.VNORMAL, - sensor.info.VSHEAR, - sensor.info.FPS, - sensor.info.FNORMAL, - sensor.info.FSHEARX, - sensor.info.FSHEARY - ) - if compute_contact: - depth, centroid, contour = sensor.read_info( - sensor.info.DEPTH, - sensor.info.CENTROID, - sensor.info.CONTOUR - ) - else: - depth = None - centroid = None - contour = None - - # Build data dictionary - data = { - 'sensor_id': sensor_id, - 'img': img if img is not None and img.size > 0 else None, - 'flow': flow if flow is not None and flow.size > 0 else None, - 'normal': normal if normal is not None and normal.size > 0 else None, - 'shear': shear if shear is not None and shear.size > 0 else None, - 'fps': fps, - 'fnormal': fnormal, - 'fshearx': fshearx, - 'fsheary': fsheary, - 'frame_count': frame_count, - } - - # Handle contact-related fields - if compute_contact: - # depth and depth_map - data['depth'] = depth if depth is not None and depth.size > 0 else None - data['depth_map'] = data['depth'] # Alias for compatibility - - # centroid - can be None or numpy array - if centroid is not None: - # Check if it's a valid array with elements - if hasattr(centroid, 'size') and centroid.size > 0: - data['centroid'] = centroid - data['centroid_found'] = True - else: - data['centroid'] = None - data['centroid_found'] = False - else: - data['centroid'] = None - data['centroid_found'] = False - - # contour - data['contour'] = contour if contour is not None and contour.size > 0 else None - else: - data['depth'] = None - data['depth_map'] = None - data['centroid'] = None - data['centroid_found'] = False - data['contour'] = None - - queue_put(data) - frame_count += 1 - - except Exception as e: - print(f"[Sensor {sensor_id}] Error: {e}") - break - - sensor.disconnect() - print(f"[Sensor {sensor_id}] Process stopped") - - -class MultiSensor: - """ - Manages multiple Sensor instances running in parallel processes. - - Each sensor runs in its own process with independent CUDA context, - enabling true parallel processing without GIL limitations. - - Example: - # Simple usage with context manager - configs = [ - {'vid_src': './s1.mp4', 'isstitch': True}, - {'vid_src': './s2.mp4', 'isstitch': True}, - ] - - with MultiSensor(configs) as multi: - for all_data in multi: - for data in all_data: - if data: - print(f"Sensor {data['sensor_id']}: FPS={data['fps']}") - cv2.imshow(f"Sensor {data['sensor_id']}", data['img']) - if cv2.waitKey(1) & 0xFF == ord('q'): - break - - # Or manual control - multi = MultiSensor(configs) - multi.start() - - data = multi.get(0) # Get from sensor 0 - all_data = multi.get_all() # Get from all sensors - - multi.stop() - """ - - def __init__(self, sensor_configs): - """ - Initialize MultiSensor manager. - - Args: - sensor_configs: List of dicts, each containing Sensor constructor args: - - vid_src: Video source (required) - int for camera ID, str for file path - - isstitch: Whether to stitch (default: True) - - cuda: Use CUDA acceleration (default: True) - - verbose: Print logs (default: False) - - config_name: Config name (default: 'prototype') - - cal_path: Calibration file path (default: None) - - Example: - configs = [ - {'vid_src': 0, 'isstitch': True, 'cuda': True}, - {'vid_src': 1, 'isstitch': True, 'cuda': True}, - ] - multi = MultiSensor(configs) - """ - import multiprocessing as mp - - self.sensor_configs = sensor_configs - self.num_sensors = len(sensor_configs) - self._mp = mp - - self.processes = [] - self.data_queues = [] - self.stop_events = [] - self._started = False - - # Create queues and events for each sensor - for i in range(self.num_sensors): - self.data_queues.append(mp.Queue(maxsize=2)) - self.stop_events.append(mp.Event()) - - def start(self): - """Start all sensor worker processes.""" - if self._started: - print("MultiSensor already started") - return - - for i, config in enumerate(self.sensor_configs): - process = self._mp.Process( - target=_sensor_worker, - args=(i, config, self.data_queues[i], self.stop_events[i]), - daemon=True - ) - process.start() - self.processes.append(process) - print(f"Started sensor {i} (PID: {process.pid})") - - self._started = True - print(f"All {len(self.processes)} sensors started") - - def get(self, sensor_idx, timeout=0.01): - """ - Get the latest data from a specific sensor. - - Args: - sensor_idx: Index of the sensor (0-based) - timeout: How long to wait for data (seconds) - - Returns: - dict: Sensor data or None if no data available - - sensor_id: int - - img: np.ndarray (H, W, 3) - - flow: np.ndarray (H, W, 2) - - normal: np.ndarray - - shear: np.ndarray - - depth: np.ndarray - - fps: float - - fnormal, fshearx, fsheary, ftor: float - - frame_count: int - """ - try: - return self.data_queues[sensor_idx].get(timeout=timeout) - except: - return None - - def get_all(self, timeout=0.01): - """ - Get the latest data from all sensors. - - Returns: - list: List of sensor data dicts (None for sensors with no new data) - """ - return [self.get(i, timeout) for i in range(self.num_sensors)] - - def __iter__(self): - """ - Iterate over sensor data. Yields all sensor data on each iteration. - - Example: - for all_data in multi: - for data in all_data: - if data: - process(data) - """ - return self - - def __next__(self): - """Get next batch of sensor data.""" - if not self._started: - raise StopIteration("MultiSensor not started") - - # Check if any process is still alive - if not any(p.is_alive() for p in self.processes): - raise StopIteration("All sensor processes have stopped") - - return self.get_all() - - def __len__(self): - """Return number of sensors.""" - return self.num_sensors - - def __getitem__(self, idx): - """Get latest data from sensor by index.""" - return self.get(idx) - - def stop(self): - """Stop all sensor processes gracefully.""" - if not self._started: - return - - print("Stopping all sensors...") - - # Signal all processes to stop - for event in self.stop_events: - event.set() - - # Wait for processes to finish - for i, process in enumerate(self.processes): - process.join(timeout=3.0) - if process.is_alive(): - print(f"Force terminating sensor {i}") - process.terminate() - - self.processes.clear() - self._started = False - print("All sensors stopped.") - - def __enter__(self): - """Context manager entry - starts sensors.""" - self.start() - return self - - def __exit__(self, exc_type, exc_val, exc_tb): - """Context manager exit - stops sensors.""" - self.stop() - return False - - def __del__(self): - """Cleanup on deletion.""" - self.stop() - diff --git a/src/orisys/StructuredGrid.py b/src/orisys/StructuredGrid.py deleted file mode 100644 index ea4d117..0000000 --- a/src/orisys/StructuredGrid.py +++ /dev/null @@ -1,408 +0,0 @@ -''' -Copyright (c) 2015, Harsh Bhatia (bhatia4@llnl.gov) -All rights reserved. - -Redistribution and use in source and binary forms, with or without -modification, are permitted provided that the following conditions are met: -1. Redistributions of source code must retain the above copyright notice, this -list of conditions and the following disclaimer. -2. Redistributions in binary form must reproduce the above copyright notice, -this list of conditions and the following disclaimer in the documentation -and/or other materials provided with the distribution. -THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND -ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED -WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE -DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE -FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL -DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER -CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, -OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -''' - -import numpy -import sys -from .timer import Timer - -# Try to import CuPy for GPU acceleration -try: - import cupy as cp - CUPY_AVAILABLE = True -except ImportError: - cp = None - CUPY_AVAILABLE = False - -# ------------------------------------------------------------------------------ -class StructuredGrid(object): - - # -------------------------------------------------------------------------- - # -------------------------------------------------------------------------- - def __init__(self, **kwargs): - - ''' - kwargs: - grid: ndarray of grid dimensions (Y,X) or (Z,Y,X) - spacings: ndarray of grid spacings (dy, dx) or (dz, dy, dx) - verbose: verbosity level - ''' - - args = kwargs.keys() - - if ('grid' not in args) or ('spacings' not in args): - raise SyntaxError("Dimensions and spacings of the grid are required") - - verbose = 0 - if 'verbose' in args: - verbose = kwargs['verbose'] - - self.dims = kwargs['grid'] - self.dx = kwargs['spacings'] - self.dim = len(self.dims) - - if self.dim != 2 and self.dim != 3: - raise ValueError("StructuredGrid works for 2D and 3D only") - - if self.dim != len(self.dx): - raise ValueError("Dimensions of spacings should match that of the grid") - - # Initialize GPU cache - self._gpu_divcurl_cache = None - - if verbose > 0: - print(' Initializing', self.dim, 'D structured grid...', end ='') - sys.stdout.flush() - mtimer = Timer() - print(' Done!', end ='') - mtimer.end() - - # -------------------------------------------------------------------------- - # -------------------------------------------------------------------------- - def divcurl(self, vfield, verbose=False): - # compute div and curl of vfield (vector field) - #if (vfield.shape[-1] != self.dim) or (vfield.shape[0:self.dim] - self.dims).any(): - if (vfield.shape[-1] != self.dim) or (vfield.shape[0:self.dim] != self.dims): - raise ValueError("Dimensions of vector field should match that of the grid") - - if verbose: - print(' Computing div and curl...', end ='') - sys.stdout.flush() - mtimer = Timer() - - if self.dim == 2: - - # self.dx = (dy,dx) - dudy, dudx = numpy.gradient(vfield[:,:,0], self.dx[0], self.dx[1]) - dvdy, dvdx = numpy.gradient(vfield[:,:,1], self.dx[0], self.dx[1]) - - numpy.add(dudx, dvdy, dudx) - numpy.subtract(dvdx, dudy, dvdx) - - if verbose: - print(' Done!', end ='') - mtimer.end() - - return (dudx, dvdx) - - elif self.dim == 3: - - # self.dx = (dz,dy,dx) - dudz, dudy, dudx = numpy.gradient(vfield[:,:,:,0], self.dx[0], self.dx[1], self.dx[2]) - dvdz, dvdy, dvdx = numpy.gradient(vfield[:,:,:,1], self.dx[0], self.dx[1], self.dx[2]) - dwdz, dwdy, dwdx = numpy.gradient(vfield[:,:,:,2], self.dx[0], self.dx[1], self.dx[2]) - - numpy.add(dudx, dvdy, dudx) - numpy.add(dudx, dvdz, dudx) - - numpy.subtract(dwdy, dvdz, dwdy) - numpy.subtract(dudz, dwdx, dudz) - numpy.subtract(dvdx, dudy, dvdx) - - if verbose: - print(' Done!', end ='') - mtimer.end() - - return (dudx, dwdy, dudz, dvdx) - - def divcurl_gpu(self, vfield, verbose=False): - """ - GPU-accelerated computation of divergence and curl of vector field. - Uses CuPy for faster gradient computation on GPU. - """ - if not CUPY_AVAILABLE: - if verbose: - print('CuPy not available, falling back to CPU') - return self.divcurl(vfield, verbose) - - if (vfield.shape[-1] != self.dim) or (vfield.shape[0:self.dim] != self.dims): - raise ValueError("Dimensions of vector field should match that of the grid") - - if verbose: - print(' Computing div and curl on GPU...', end ='') - sys.stdout.flush() - mtimer = Timer() - - if self.dim == 2: - # Transfer to GPU - vfield_gpu = cp.asarray(vfield, dtype=cp.float32) - - # Compute gradients on GPU - dudy, dudx = cp.gradient(vfield_gpu[:,:,0], self.dx[0], self.dx[1]) - dvdy, dvdx = cp.gradient(vfield_gpu[:,:,1], self.dx[0], self.dx[1]) - - # Compute div and curl on GPU - div_gpu = dudx + dvdy - curl_gpu = dvdx - dudy - - # Transfer back to CPU - div = cp.asnumpy(div_gpu) - curl = cp.asnumpy(curl_gpu) - - if verbose: - print(' Done!', end ='') - mtimer.end() - - return (div, curl) - - elif self.dim == 3: - # Transfer to GPU - vfield_gpu = cp.asarray(vfield, dtype=cp.float32) - - # Compute gradients on GPU - dudz, dudy, dudx = cp.gradient(vfield_gpu[:,:,:,0], self.dx[0], self.dx[1], self.dx[2]) - dvdz, dvdy, dvdx = cp.gradient(vfield_gpu[:,:,:,1], self.dx[0], self.dx[1], self.dx[2]) - dwdz, dwdy, dwdx = cp.gradient(vfield_gpu[:,:,:,2], self.dx[0], self.dx[1], self.dx[2]) - - # Compute div and curl components on GPU - div_gpu = dudx + dvdy + dwdz - curlu_gpu = dwdy - dvdz - curlv_gpu = dudz - dwdx - curlw_gpu = dvdx - dudy - - # Transfer back to CPU - div = cp.asnumpy(div_gpu) - curlu = cp.asnumpy(curlu_gpu) - curlv = cp.asnumpy(curlv_gpu) - curlw = cp.asnumpy(curlw_gpu) - - if verbose: - print(' Done!', end ='') - mtimer.end() - - return (div, curlu, curlv, curlw) - - def curl3D(self, vfield, verbose=False): - - #if (vfield.shape[-1] != self.dim) or (vfield.shape[0:self.dim] - self.dims).any(): - if (vfield.shape[-1] != self.dim) or (vfield.shape[0:self.dim] != self.dims): - raise ValueError("Dimensions of vector field should match that of the grid") - - if self.dim != 3: - raise ValueError("curl3D works only for 2D") - - if verbose: - print(' Computing curl...', end ='') - sys.stdout.flush() - mtimer = Timer() - - # self.dx = (dz,dy,dx) - dudz, dudy, dudx = numpy.gradient(vfield[:,:,:,0], self.dx[0], self.dx[1], self.dx[2]) - dvdz, dvdy, dvdx = numpy.gradient(vfield[:,:,:,1], self.dx[0], self.dx[1], self.dx[2]) - dwdz, dwdy, dwdx = numpy.gradient(vfield[:,:,:,2], self.dx[0], self.dx[1], self.dx[2]) - - numpy.subtract(dwdy, dvdz, dwdy) - numpy.subtract(dudz, dwdx, dudz) - numpy.subtract(dvdx, dudy, dvdx) - - if verbose: - print(' Done!', end ='') - mtimer.end() - - return (dwdy, dudz, dvdx) - - def curl3D_gpu(self, vfield, verbose=False): - """ - GPU-accelerated computation of curl for 3D vector field. - Uses CuPy for faster gradient computation on GPU. - """ - if not CUPY_AVAILABLE: - if verbose: - print('CuPy not available, falling back to CPU') - return self.curl3D(vfield, verbose) - - if (vfield.shape[-1] != self.dim) or (vfield.shape[0:self.dim] != self.dims): - raise ValueError("Dimensions of vector field should match that of the grid") - - if self.dim != 3: - raise ValueError("curl3D works only for 3D") - - if verbose: - print(' Computing curl on GPU...', end ='') - sys.stdout.flush() - mtimer = Timer() - - # Transfer to GPU - vfield_gpu = cp.asarray(vfield, dtype=cp.float32) - - # Compute gradients on GPU - dudz, dudy, dudx = cp.gradient(vfield_gpu[:,:,:,0], self.dx[0], self.dx[1], self.dx[2]) - dvdz, dvdy, dvdx = cp.gradient(vfield_gpu[:,:,:,1], self.dx[0], self.dx[1], self.dx[2]) - dwdz, dwdy, dwdx = cp.gradient(vfield_gpu[:,:,:,2], self.dx[0], self.dx[1], self.dx[2]) - - # Compute curl components on GPU - curl_u = dwdy - dvdz - curl_v = dudz - dwdx - curl_w = dvdx - dudy - - # Transfer back to CPU - curlu = cp.asnumpy(curl_u) - curlv = cp.asnumpy(curl_v) - curlw = cp.asnumpy(curl_w) - - if verbose: - print(' Done!', end ='') - mtimer.end() - - return (curlu, curlv, curlw) - - def rotated_gradient(self, sfield, verbose=False): - - if (sfield.shape != self.dims): - #if (sfield.shape - self.dims).any(): - raise ValueError("Dimensions of scalar field should match that of the grid") - - if self.dim != 2: - raise ValueError("rotated_gradient works only for 2D") - - if verbose: - print(' Computing rotated gradient...', end ='') - sys.stdout.flush() - mtimer = Timer() - - ddy, ddx = numpy.gradient(sfield, self.dx[0], self.dx[1]) - ddy *= -1.0 - - grad = numpy.stack((ddy, ddx), axis=-1) - - if verbose: - print(' Done!', end ='') - mtimer.end() - - return grad - - def rotated_gradient_gpu(self, sfield, verbose=False): - """ - GPU-accelerated computation of rotated gradient for 2D scalar field. - Uses CuPy for faster gradient computation on GPU. - """ - if not CUPY_AVAILABLE: - if verbose: - print('CuPy not available, falling back to CPU') - return self.rotated_gradient(sfield, verbose) - - if (sfield.shape != self.dims): - raise ValueError("Dimensions of scalar field should match that of the grid") - - if self.dim != 2: - raise ValueError("rotated_gradient works only for 2D") - - if verbose: - print(' Computing rotated gradient on GPU...', end ='') - sys.stdout.flush() - mtimer = Timer() - - # Transfer to GPU - sfield_gpu = cp.asarray(sfield, dtype=cp.float32) - - # Compute gradients on GPU - ddy, ddx = cp.gradient(sfield_gpu, self.dx[0], self.dx[1]) - ddy *= -1.0 - - # Stack and transfer back to CPU - grad_gpu = cp.stack((ddy, ddx), axis=-1) - grad = cp.asnumpy(grad_gpu) - - if verbose: - print(' Done!', end ='') - mtimer.end() - - return grad - - def gradient(self, sfield, verbose=False): - - if (sfield.shape != self.dims): - #if (sfield.shape - self.dims).any(): - raise ValueError("Dimensions of scalar field should match that of the grid") - - if verbose: - print(' Computing gradient...', end ='') - sys.stdout.flush() - mtimer = Timer() - - if self.dim == 2: - - # self.dx = (dy,dx) - ddy, ddx = numpy.gradient(sfield, self.dx[0], self.dx[1]) - grad = numpy.stack((ddx, ddy), axis = -1) - - elif self.dim == 3: - - # self.dx = (dz,dy,dx) - ddz, ddy, ddx = numpy.gradient(sfield, self.dx[0], self.dx[1], self.dx[2]) - grad = numpy.stack((ddx, ddy, ddz), axis = -1) - - if verbose: - print(' Done!', end ='') - mtimer.end() - - return grad - - def gradient_gpu(self, sfield, verbose=False): - """ - GPU-accelerated computation of gradient for scalar field. - Uses CuPy for faster gradient computation on GPU. - """ - if not CUPY_AVAILABLE: - if verbose: - print('CuPy not available, falling back to CPU') - return self.gradient(sfield, verbose) - - if (sfield.shape != self.dims): - raise ValueError("Dimensions of scalar field should match that of the grid") - - if verbose: - print(' Computing gradient on GPU...', end ='') - sys.stdout.flush() - mtimer = Timer() - - # Transfer to GPU - sfield_gpu = cp.asarray(sfield, dtype=cp.float32) - - if self.dim == 2: - # Compute gradients on GPU - ddy, ddx = cp.gradient(sfield_gpu, self.dx[0], self.dx[1]) - grad_gpu = cp.stack((ddx, ddy), axis=-1) - - elif self.dim == 3: - # Compute gradients on GPU - ddz, ddy, ddx = cp.gradient(sfield_gpu, self.dx[0], self.dx[1], self.dx[2]) - grad_gpu = cp.stack((ddx, ddy, ddz), axis=-1) - - # Transfer back to CPU - grad = cp.asnumpy(grad_gpu) - - if verbose: - print(' Done!', end ='') - mtimer.end() - - return grad - - def clear_gpu_cache(self): - """Clear GPU cache to free GPU memory""" - if self._gpu_divcurl_cache is not None: - self._gpu_divcurl_cache = None - if CUPY_AVAILABLE: - cp.get_default_memory_pool().free_all_blocks() - - # -------------------------------------------------------------------------- - # -------------------------------------------------------------------------- diff --git a/src/orisys/UnstructuredGrid.py b/src/orisys/UnstructuredGrid.py deleted file mode 100644 index 698247a..0000000 --- a/src/orisys/UnstructuredGrid.py +++ /dev/null @@ -1,424 +0,0 @@ -''' -Copyright (c) 2015, Harsh Bhatia (bhatia4@llnl.gov) -All rights reserved. - -Redistribution and use in source and binary forms, with or without -modification, are permitted provided that the following conditions are met: -1. Redistributions of source code must retain the above copyright notice, this -list of conditions and the following disclaimer. -2. Redistributions in binary form must reproduce the above copyright notice, -this list of conditions and the following disclaimer in the documentation -and/or other materials provided with the distribution. -THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND -ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED -WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE -DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE -FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL -DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER -CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, -OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -''' - -import sys -import numpy -from scipy import spatial -from .timer import Timer - -class UnstructuredGrid(object): - - '''Class to support nHHD on unstructured grids (triangular and tetrahedral)''' - - # -------------------------------------------------------------------------- - # -------------------------------------------------------------------------- - def need_volumes(self, verbose=False): - - ''' Compute volumes/areas for vertices, simplices, and corners''' - - if self.pvolumes.shape != (0,): - return self.pvolumes - - def tri_area(o,a,b): - return numpy.linalg.norm(numpy.cross(b-o,a-o)) / 2.0 - - def tet_volume(o,a,b,c): - return numpy.abs(numpy.dot(a-o,numpy.cross(b-o,c-o))) / 6.0 - - if verbose: - print(' Computing point areas/volumes...', end ='') - sys.stdout.flush() - mtimer = Timer() - - self.svolumes = numpy.zeros(self.nsimplices,) - self.pvolumes = numpy.zeros(self.nvertices,) - self.cvolumes = numpy.zeros((self.nsimplices,self.dim+1)) - - # triangulation - if self.dim == 2: - - for sidx in xrange(self.nsimplices): - - simp = self.simplices[sidx] - verts = self.vertices[simp] - - self.svolumes[sidx] = tri_area(verts[0], verts[1], verts[2]) - - # edge weights - e = [ verts[2]-verts[1], verts[0]-verts[2], verts[1]-verts[0] ] - l2 = [ numpy.dot(e[0],e[0]), numpy.dot(e[1],e[1]), numpy.dot(e[2],e[2]) ] - ew = [ l2[0]*(l2[1]+l2[2]-l2[0]), l2[1]*(l2[2]+l2[0]-l2[1]),l2[2]*(l2[0]+l2[1]-l2[2]) ] - - # corner areas - if (ew[0] <= 0): - self.cvolumes[sidx,1] = -0.25 * l2[2] * self.svolumes[sidx] / numpy.dot(e[0], e[2]) - self.cvolumes[sidx,2] = -0.25 * l2[1] * self.svolumes[sidx] / numpy.dot(e[0], e[1]) - self.cvolumes[sidx,0] = self.svolumes[sidx] - self.cvolumes[sidx,1] - self.cvolumes[sidx,2] - elif (ew[1] <= 0): - self.cvolumes[sidx,2] = -0.25 * l2[0] * self.svolumes[sidx] / numpy.dot(e[1], e[0]); - self.cvolumes[sidx,0] = -0.25 * l2[2] * self.svolumes[sidx] / numpy.dot(e[1], e[2]); - self.cvolumes[sidx,1] = self.svolumes[sidx] - self.cvolumes[sidx,2] - self.cvolumes[sidx,0]; - elif (ew[2] <= 0): - self.cvolumes[sidx,0] = -0.25 * l2[1] * self.svolumes[sidx] / numpy.dot(e[2], e[1]); - self.cvolumes[sidx,1] = -0.25 * l2[0] * self.svolumes[sidx] / numpy.dot(e[2], e[0]); - self.cvolumes[sidx,2] = self.svolumes[sidx] - self.cvolumes[sidx,0] - self.cvolumes[sidx,1]; - else: - ewscale = 0.5 * self.svolumes[sidx] / (ew[0] + ew[1] + ew[2]) - for d in xrange(3): - self.cvolumes[sidx,d] = ewscale * (ew[(d+1)%3] + ew[(d+2)%3]) - - self.pvolumes[simp[0]] += self.cvolumes[sidx,0] - self.pvolumes[simp[1]] += self.cvolumes[sidx,1] - self.pvolumes[simp[2]] += self.cvolumes[sidx,2] - - # tetrahedralization - elif self.sdim == 3: - - raise ValueError('TODO: pvolumes for 3D') - for sidx in xrange(self.nsimplices): - - simp = self.simplices[sidx] - verts = self.vertices[simp] - - self.svolumes[sidx] = tet_volume(verts[0], verts[1], verts[2], verts[3]) - - for v in simp: - self.pvolumes[v] += self.svolumes[sidx] / 4.0 - - - if verbose: - print(' Done!', end ='') - mtimer.end() - - return self.pvolumes - - def need_adjacentfaces(self, verbose=False): - - ''' - Find adjacent faces for each vertex - as list of lists - ''' - - if len(self.adjacent_faces) != 0: - return self.adjacent_faces - - if verbose: - print(' Computing adjacent_faces...', end ='') - sys.stdout.flush() - mtimer = Timer() - - numadjacentfaces = numpy.zeros(self.nvertices, dtype=int) - for f in self.simplices: - for i in xrange(3): - numadjacentfaces[f[i]] += 1 - - # can be optimized further by avoiding "append"? - self.adjacent_faces = [[] for _ in xrange(self.nvertices)] - for fidx in xrange(self.nsimplices): - for i in xrange(3): - self.adjacent_faces[self.simplices[fidx, i]].append(fidx) - - if verbose: - print('Done!', end ='') - mtimer.end() - - return self.adjacent_faces - - def need_acrossedge(self, verbose=False): - - ''' - Find adjacent faces for each face (across each edge) - as ndarray of ints: shape (nsimplex, 3) - -1 denotes a face on the boundary (no face across edge) - ''' - if self.across_edge.shape != (0,0): - return self.across_edge - - self.need_adjacentfaces(verbose) - - if verbose: - print(' Computing across_edge...', end ='') - sys.stdout.flush() - mtimer = Timer() - - self.across_edge = -1 * numpy.ones((self.nsimplices, 3), dtype=numpy.int) - - for fidx in xrange(self.nsimplices): - for i in xrange(3): - - if self.across_edge[fidx, i] != -1: - continue - - v1 = self.simplices[fidx, (i+1)%3] - v2 = self.simplices[fidx, (i+2)%3] - - for other in self.adjacent_faces[v1]: - - if other == fidx: - continue - - if other not in self.adjacent_faces[v2]: - continue - - oface = self.simplices[other] - j = numpy.where(oface == v1)[0] - j = (j+1)%3 - - if oface[(j+1)%3] != v2: - continue - - self.across_edge[fidx, i] = other - self.across_edge[other, j] = fidx - - if verbose: - print(' Done!', end ='') - mtimer.end() - return self.across_edge - - def need_boundary(self, verbose=False): - - ''' - Find boundary of the triangulation - as collection of boundary edeges as an array of [face, k] - ''' - - if self.bedges.shape != (0,0): - return self.bedges - - self.need_acrossedge(verbose) - - if verbose: - print(' Computing the boundary...', end ='') - sys.stdout.flush() - mtimer = Timer() - - # find all boundary faces and edges - bfaces = [fidx for fidx in xrange(self.nsimplices) if -1 in self.across_edge[fidx]] - bedges = [] - - for fidx in bfaces: - face = self.simplices[fidx] - nbrs = self.across_edge[fidx] - bedges.extend( [[fidx, k] for k in xrange(3) if nbrs[k] == -1] ) - - self.bedges = numpy.array(bedges) - if verbose: - print(' Done! found', self.bedges.shape[0], 'boundary edges', end ='') - mtimer.end() - - return self.bedges - - - # -------------------------------------------------------------------------- - # -------------------------------------------------------------------------- - def __init__(self, **kwargs): - - ''' - kwargs: - vertices: ndarray of shape (nverts, dim) # dim = 2,3 - simplices: ndarray of shape (nfaces, dim+1) - verbose: verbosity level - ''' - - args = kwargs.keys() - - if 'vertices' not in args: - raise SyntaxError("Mesh object needs vertex data") - - verbose = 2 - if 'verbose' in args: - verbose = kwargs['verbose'] - - self.vertices = kwargs['vertices'] - self.dim = self.vertices.shape[1] - self.nvertices = self.vertices.shape[0] - - if self.dim != 2 and self.dim != 3: - raise SyntaxError("Mesh object works for 2D and 3D only") - - if verbose > 0: - print(' Initializing', self.dim, 'D mesh with', self.nvertices, 'vertices...', end ='') - sys.stdout.flush() - mtimer = Timer() - - if verbose > 1: - print('') - - # create simplices if needed - if 'simplices' in args: - self.Delaunay = None - self.simplices = kwargs['simplices'] - if verbose > 1: - print(' got', self.simplices.shape[0], 'simplices') - else: - if verbose > 1: - print(' creating Delaunay mesh...') - sys.stdout.flush() - - self.Delaunay = spatial.Delaunay(self.vertices) - self.simplices = self.Delaunay.simplices - - if verbose > 1: - print(' Done! created', self.simplices.shape[0], 'simplices') - sys.stdout.flush() - - self.nsimplices = self.simplices.shape[0] - if self.dim != self.simplices.shape[1]-1: - raise SyntaxError("Dimension mismatch! pdim = "+str(pdim)+" and sdim = "+str(sdim)+" do not match!") - - self.adjacent_faces = [] - self.across_edge = numpy.empty((0,0)) - self.bedges = numpy.empty((0,0)) - self.pvolumes = numpy.empty(0) - - self.need_volumes(verbose > 1) - self.need_boundary(verbose > 1) - #self.need_meshmatrices(verbose > 1) - - if verbose > 0: - print(' Done!', end ='') - mtimer.end() - - # -------------------------------------------------------------------------- - # -------------------------------------------------------------------------- - def divcurl(self, vfield, verbose=False): - - if vfield.shape != (self.nsimplices, self.dim): - print('UnstructuredGrid.divcurl', vfield.shape, self.dim, self.vertices.shape, self.simplices.shape) - raise ValueError("UnstructuredGrid requires a valid-dimensional vector field") - - if verbose: - print(' Computing divcurl...', end ='') - sys.stdout.flush() - mtimer = Timer() - - div = numpy.zeros(self.nvertices) - curlw = numpy.zeros(self.nvertices) - - # for each face - for sidx in xrange(self.nsimplices): - - simp = self.simplices[sidx] - - for k in xrange(3): - - v = simp[k] - a = simp[ (k+1)%3 ] - b = simp[ (k+2)%3 ] - - # normal and tangent vectors - tvec = self.vertices[b] - self.vertices[a] # counterclockwise - nvec = numpy.array([-tvec[1], tvec[0]]) # inward - - dn = numpy.dot(nvec, vfield[sidx]) - tn = numpy.dot(tvec, vfield[sidx]) - - div[v] += dn - curlw[v] += tn - - # fix for boundary edges - for bedge in self.bedges: - - sidx = bedge[0] - eidx = bedge[1] - - a = self.simplices[sidx][(eidx+1)%3] - b = self.simplices[sidx][(eidx+2)%3] - - tvec = self.vertices[b] - self.vertices[a] # counterclockwise - nvec = numpy.array([-tvec[1], tvec[0]]) # inward - - dn = numpy.dot(nvec, vfield[sidx]) - dt = numpy.dot(tvec, vfield[sidx]) - - div[a] += dn - div[b] += dn - curlw[a] += tn - curlw[b] += tn - - div = -0.5 * div - curlw = 0.5 * curlw - - if verbose: - print(' Done!', end ='') - mtimer.end() - #print('\tdiv =', div.shape, div.min(), div.max()) - #print('\tcurlw =', curlw.shape, curlw.min(), curlw.max()) - - return (div, curlw) - - def gradient(self, sfield, verbose=False): - - if sfield.shape[0] != self.nvertices: - print('UnstructuredGrid.gradient', self.dim, sfield.shape, self.vertices.shape, self.simplices.shape) - raise ValueError("UnstructuredGrid requires a valid-dimensional vector field") - - if verbose: - print(' Computing gradient...', end ='') - sys.stdout.flush() - mtimer = Timer() - - grad = numpy.zeros((self.nsimplices, self.dim)) - - # for 2D - for sidx in xrange(self.nsimplices): - - deb = sidx==0 - - simp = self.simplices[sidx] - f = 0.5 / self.svolumes[sidx] - - if deb: - print(sidx, simp, self.svolumes[sidx], f) - - for k in xrange(3): - - v = simp[k] - a = simp[ (k+1)%3 ] - b = simp[ (k+2)%3 ] - - # normal and tangent vectors - tvec = self.vertices[b] - self.vertices[a] # counterclockwise - nvec = numpy.array([-tvec[1], tvec[0]]) # inward - - grad[sidx] += f * sfield[v] * nvec - - if deb: - print('\t', k, nvec, sfield[v],grad[sidx]) - - if verbose: - print(' Done!', end ='') - mtimer.end() - return grad - - def rotated_gradient(self, sfield): - - rgrad = self.gradient(sfield) - - rgrad[:,[0, 1]] = rgrad[:,[1, 0]] - rgrad[:,0] *= -1.0 - - return rgrad - - # -------------------------------------------------------------------------- - # -------------------------------------------------------------------------- diff --git a/src/orisys/__init__.py b/src/orisys/__init__.py index edfbfca..a3a6ffb 100644 --- a/src/orisys/__init__.py +++ b/src/orisys/__init__.py @@ -1,5 +1,5 @@ """ -Orisys - Optical Tactile Sensor SDK +orisys - Optical Tactile Sensor SDK A Python library for processing and analyzing optical tactile sensor data, including image stitching, optical flow tracking, and deformation decomposition. @@ -9,7 +9,21 @@ __version__ = "0.2.0" __author__ = "Orisys" from .sensor import Sensor -from .MultiSensor import MultiSensor from . import util +from . import configs -__all__ = ['Sensor', 'MultiSensor', 'util'] +# Re-export list_configs as list_config for convenience +list_config = configs.list_config +export_config = configs.export_config + +__all__ = [ + # class + 'Sensor', + + # util functions + 'util', + + # configs utilities functions + 'list_config', + 'export_config' +] diff --git a/src/orisys/configs/__init__.py b/src/orisys/configs/__init__.py index b4f5381..97d8b65 100644 --- a/src/orisys/configs/__init__.py +++ b/src/orisys/configs/__init__.py @@ -1,7 +1,57 @@ import json +import sys from pathlib import Path -CONFIG_DIR = Path(__file__).parent +def _get_config_dir(): + """Resolve config directory, robustly handling Nuitka compiled paths.""" + + # 1. Start with the directory of this file + current_path = Path(__file__).parent + + # 2. Check strict relative location (standard case) + # If we are in src/orisys/configs/__init__.py, configs is current_path + if (current_path / "prototype.json").exists(): + return current_path + + # 3. Walk up the directory tree to find a 'configs' folder containing 'prototype.json' + # This handles nested virtual paths like .../orisys/orisys/configs + # We look for a physical directory that matches our expectation + + search_path = current_path + # Limit depth to avoid infinite loops, though unlikely with .parent + for _ in range(5): + # Check if 'configs' exists in this directory + candidate = search_path / "configs" + if candidate.exists() and candidate.is_dir(): + if (candidate / "prototype.json").exists(): + return candidate + + # Also check if the current directory ITSELF is the configs dir we missed + # (e.g. if we are somehow inside configs but missed check 1) + if search_path.name == 'configs' and (search_path / "prototype.json").exists(): + return search_path + + if search_path.parent == search_path: # Reached root + break + search_path = search_path.parent + + # 4. Try finding 'orisys' package root via sys.modules (if loaded) + try: + if 'orisys' in sys.modules: + mod = sys.modules['orisys'] + if hasattr(mod, '__file__') and mod.__file__: + pkg_root = Path(mod.__file__).parent + # pkg_root is .../site-packages/orisys + candidate = pkg_root / "configs" + if candidate.exists() and (candidate / "prototype.json").exists(): + return candidate + except Exception: + pass + + # Fallback to current path if nothing else works + return current_path + +CONFIG_DIR = _get_config_dir() def load_config(config_name_or_path): """ @@ -56,7 +106,7 @@ def load_config(config_name_or_path): with open(package_config_path, 'r', encoding='utf-8') as f: return json.load(f) -def list_configs(): +def list_config(): """ List all available JSON configuration files in the configs directory. @@ -67,6 +117,29 @@ def list_configs(): config_files = sorted(CONFIG_DIR.glob("*.json")) return [f.stem for f in config_files] +def export_config(config_name = None, path="./config/exported_config.json"): + """ + Export configuration to JSON file with readable formatting. + Arrays/lists are kept on a single line, while nested objects are indented. + + Args: + config: Configuration dictionary to export + path: Output file path (default: "./config/exported_config.json") + + Returns: + The exported config dictionary + """ + config = load_config(config_name) + Path(path).parent.mkdir(parents=True, exist_ok=True) + + with open(path, 'w', encoding='utf-8') as f: + formatted = _format_json_value(config) + f.write(formatted) + f.write('\n') # Add trailing newline + + print(f"Config name: {config_name} exported to {path}") + return config + def _format_json_value(value, indent_level=0, indent_str=" "): """ Recursively format JSON value with arrays/lists on single line. @@ -96,26 +169,4 @@ def _format_json_value(value, indent_level=0, indent_str=" "): return "[" + ", ".join(items) + "]" else: # Primitive types (string, number, bool, None) - return json.dumps(value, ensure_ascii=False) - -def export_config(config, path="./config/exported_config.json"): - """ - Export configuration to JSON file with readable formatting. - Arrays/lists are kept on a single line, while nested objects are indented. - - Args: - config: Configuration dictionary to export - path: Output file path (default: "./config/exported_config.json") - - Returns: - The exported config dictionary - """ - Path(path).parent.mkdir(parents=True, exist_ok=True) - - with open(path, 'w', encoding='utf-8') as f: - formatted = _format_json_value(config) - f.write(formatted) - f.write('\n') # Add trailing newline - - return config - \ No newline at end of file + return json.dumps(value, ensure_ascii=False) \ No newline at end of file diff --git a/src/orisys/getDensity.py b/src/orisys/getDensity.py deleted file mode 100644 index c0730b4..0000000 --- a/src/orisys/getDensity.py +++ /dev/null @@ -1,356 +0,0 @@ -""" -Density estimation module for optical flow fields. - -This module computes density maps from optical flow vectors using Gaussian kernel -density estimation. It supports both CPU (NumPy) and CUDA (CuPy) implementations. -""" -import cupy as cp -import cv2 -import numpy as np - -# Global variable to store the baseline/zero level of density for invalid area correction -zero_level = None - -def getFilteredDensity(flow, use_cuda): - """ - Compute filtered density map from optical flow field. - - This function computes a density map using Gaussian kernel density estimation - and then applies filtering to correct invalid areas. - - Args: - flow: Optical flow field with shape (height, width, 2) containing flow vectors - use_cuda: If True, use CUDA-accelerated computation (requires CuPy and CUDA) - If False, use CPU-based NumPy implementation - - Returns: - np.ndarray: Filtered density map with shape (height, width), where negative - values indicate regions of flow divergence - """ - # Gaussian kernel parameters for density estimation - sigma = 3.0 # Standard deviation of Gaussian kernel - r_max = 9.0 # Maximum radius for Gaussian contribution (cutoff distance) - - # Compute density map using either CUDA or CPU implementation - if use_cuda: - density = getDensityCuda(flow, sigma, r_max) - else: - density = getDensity(flow, sigma, r_max) - - # Apply filtering to correct invalid areas (regions with unreliable flow) - density = correctInvalidArea(density) - - return density - -def correctInvalidArea(density): - """ - Correct invalid areas in the density map. - - Invalid areas are regions where the density falls below a threshold (zero_level), - which typically indicates unreliable or missing flow data. This function: - 1. Uses adaptive zero_level that updates gradually to handle changing conditions - 2. Identifies invalid regions using adaptive thresholding (below zero_level - margin) - 3. Fills invalid areas using morphological closing for smoother transitions - 4. Applies stronger edge-preserving smoothing to eliminate discontinuities - - Args: - density: Input density map with shape (height, width) - - Returns: - np.ndarray: Corrected density map with invalid areas filled and smoothed, - with improved continuity to reduce black dots/discontinuities - """ - global zero_level - - # Compute adaptive zero level from current frame's central region - # This adapts to changing density values when sensor is pressed - current_zero_level = np.mean(density[60:-60, 60:-60]) - - # Initialize or update zero_level with exponential moving average for stability - if zero_level is None: - zero_level = current_zero_level - else: - # Use exponential moving average: 90% old, 10% new (adapts slowly) - zero_level = 0.9 * zero_level + 0.1 * current_zero_level - - # Use a threshold slightly below zero level to avoid over-correction - # This prevents valid but slightly lower density regions from being marked invalid - threshold = zero_level - np.std(density[60:-60, 60:-60]) * 0.5 - invalid_ind = twoClustering(density, threshold) - - # Fill invalid areas using morphological closing to smooth transitions - # This fills small holes and provides better continuity - if np.any(invalid_ind > 0): - # Create a temporary copy with invalid areas filled with zero_level - density_filled = density.copy() - density_filled[np.where(invalid_ind > 0)] = zero_level - - # Apply morphological closing to fill small holes and smooth boundaries - # This helps connect nearby valid regions and eliminate black dots - kernel = np.ones((5, 5), np.float32) - density_closed = cv2.morphologyEx(density_filled, cv2.MORPH_CLOSE, kernel, iterations=1) - - # Use closed version only for invalid regions, preserve original valid values - density[np.where(invalid_ind > 0)] = density_closed[np.where(invalid_ind > 0)] - - # Apply edge-preserving smoothing using guided filter with stronger smoothing - # Increased radius and epsilon for better continuity and to eliminate discontinuities - try: - density = cv2.ximgproc.guidedFilter(density, density, 15, 0.01) - except AttributeError: - # Fallback to stronger Gaussian blur if ximgproc is not available - # Use larger kernel for better smoothing of discontinuities - density = cv2.GaussianBlur(density, (21, 21), 0) - - return density - - -def twoClustering(density, threshold): - """ - Perform binary clustering/segmentation of density map based on threshold. - - This function creates a binary mask where pixels below the threshold are - marked as invalid (value 1) and pixels above threshold are valid (value 0). - - Args: - density: Input density map - threshold: Threshold value for segmentation - - Returns: - np.ndarray: Binary mask with dtype uint8, where 1 indicates invalid areas - (density < threshold) and 0 indicates valid areas - """ - cluster = np.zeros_like(density, dtype='uint8') - cluster[np.where(density < threshold)] = 1 - return cluster - - - -def getDensity(flow, sigma, r_max): - """ - Compute Gaussian density map from optical flow vectors. - Pure NumPy implementation without freud dependency. - Vectorized for better performance. - - Args: - flow: Optical flow field with shape (height, width, 2) - sigma: Gaussian kernel standard deviation - r_max: Maximum radius for Gaussian contribution - - Returns: - Negative density map with same shape as flow field (height, width) - """ - # Pre-compute constants - two_sigma2 = 2.0 * sigma * sigma - gaussian_norm = 1.0 / (two_sigma2 * np.pi) - - # Get flow dimensions - height, width = flow.shape[:2] - flow_x = flow[:, :, 0].astype(np.float32) - flow_y = flow[:, :, 1].astype(np.float32) - - # Create coordinate grids - x_coords = np.arange(width, dtype=np.float32) - y_coords = np.arange(height, dtype=np.float32) - grid_x, grid_y = np.meshgrid(x_coords, y_coords) - - # Compute destination positions after flow displacement - dest_x = grid_x + flow_x - dest_y = grid_y + flow_y - - # Initialize result array - result = np.zeros((height, width), dtype=np.float32) - - # Compute density by iterating over each source pixel - # For each pixel, apply its flow vector to get destination position, - # then accumulate Gaussian contributions to nearby grid points - for y in range(height): - for x in range(width): - # Get the destination position after applying flow displacement - dest_x_val = dest_x[y, x] - dest_y_val = dest_y[y, x] - - # Determine the bounding box of grid points affected by this flow vector - # These are grid points within r_max radius of the destination - min_x = max(0, int(np.floor(dest_x_val - r_max))) - max_x = min(width, int(np.ceil(dest_x_val + r_max)) + 1) - min_y = max(0, int(np.floor(dest_y_val - r_max))) - max_y = min(height, int(np.ceil(dest_y_val + r_max)) + 1) - - # Skip if bounding box is empty or out of bounds - if min_x >= max_x or min_y >= max_y: - continue - - # Create local coordinate grid for the affected region - local_x = np.arange(min_x, max_x, dtype=np.float32) - local_y = np.arange(min_y, max_y, dtype=np.float32) - local_grid_x, local_grid_y = np.meshgrid(local_x, local_y) - - # Compute squared distances from destination to each grid point - dx = local_grid_x - dest_x_val - dy = local_grid_y - dest_y_val - distance_sq = dx * dx + dy * dy - - # Create mask for points within r_max radius (circular neighborhood) - mask = distance_sq <= (r_max * r_max) - - # Compute Gaussian kernel values: exp(-d²/(2σ²)) / (2πσ²) - # Higher values for closer points, decreasing with distance - gaussian_vals = np.exp(-distance_sq / two_sigma2) * gaussian_norm - gaussian_vals[~mask] = 0 # Zero out contributions beyond r_max - - # Accumulate negative density (negative indicates flow divergence) - # Each flow vector contributes to density at its destination region - result[min_y:max_y, min_x:max_x] -= gaussian_vals - - return result - - -class CudaArrayInterface: - """ - Interface class to convert OpenCV CUDA GpuMat to CuPy-compatible array interface. - - This allows seamless conversion between OpenCV's CUDA memory objects and CuPy arrays - without copying data, by implementing the CUDA array interface protocol. - """ - def __init__(self, gpu_mat): - """ - Initialize CUDA array interface from OpenCV GpuMat. - - Args: - gpu_mat: OpenCV CUDA GpuMat object - """ - w, h = gpu_mat.size() - - # Map OpenCV data types to NumPy/CuPy type strings - type_map = { - cv2.CV_8U: "u1", # unsigned 8-bit - cv2.CV_8S: "i1", # signed 8-bit - cv2.CV_16U: "u2", # unsigned 16-bit - cv2.CV_16S: "i2", # signed 16-bit - cv2.CV_32S: "i4", # signed 32-bit integer - cv2.CV_32F: "f4", # 32-bit float - cv2.CV_64F: "f8", # 64-bit float - } - - # Implement CUDA array interface protocol for zero-copy conversion - self.__cuda_array_interface__ = { - "version": 2, - "shape": (h, w), - "data": (gpu_mat.cudaPtr(), False), # CUDA memory pointer, not writable - "typestr": type_map[gpu_mat.type()], - "strides": (gpu_mat.step, gpu_mat.elemSize()), # Row stride and element size - } - - -# CUDA kernel code for accelerated density computation -# This kernel computes Gaussian density contributions in parallel on GPU -func_string = r""" - extern "C"{ - // Compute 2D Gaussian kernel value - // Args: - // twoSigma2: 2 * sigma^2 (precomputed for efficiency) - // dx, dy: Distance components from center - // Returns: Gaussian kernel value normalized by (2*pi*sigma^2) - __device__ float gaussian(float twoSigma2, float dx, float dy) - { - float pi = 3.1415926535897932384626; - float distance = dx*dx+dy*dy; // Squared distance - float result = expf(-distance/twoSigma2); // Gaussian exponential term - result = result/(twoSigma2*pi); // Normalization factor - return result; - } - - // CUDA kernel to compute density map from optical flow - // Each thread processes one pixel in the flow field - // Args: - // flowX, flowY: Input flow field components (flattened) - // result: Output density map (flattened, initialized to zero) - // x_size, y_size: Dimensions of the flow field - // twoSigma2: Precomputed 2 * sigma^2 - // r_max: Maximum radius for Gaussian contribution - __global__ void density(const float *flowX, const float *flowY, float *result, - int x_size, int y_size, float twoSigma2, float r_max) - { - // Get thread coordinates (each thread handles one source pixel) - int x = threadIdx.x; // Column index - int y = blockIdx.x; // Row index - - // Compute destination position after applying flow displacement - float dest_x = float(x) + flowX[y*x_size+x]; - float dest_y = float(y) + flowY[y*x_size+x]; - - // Iterate over nearby grid points within r_max radius - for(int ind_x = ceilf(dest_x - r_max); ind_x < ceilf(dest_x + r_max); ++ind_x) - for(int ind_y = ceilf(dest_y - r_max); ind_y < ceilf(dest_y + r_max); ++ind_y) - { - // Bounds checking - if(ind_x >= x_size || ind_x < 0 || ind_y >= y_size || ind_y < 0) - { - continue; - } - - // Compute distance from destination to grid point - float dx = float(ind_x) - dest_x; - float dy = float(ind_y) - dest_y; - - // Accumulate Gaussian contribution (negative for divergence representation) - atomicAdd(result + ind_y*x_size + ind_x, -gaussian(twoSigma2, dx, dy)); - } - } - } - """ -# Compile CUDA kernel into CuPy RawModule for execution -funcupy = cp.RawModule(code=func_string).get_function("density") - - -def getDensityCuda(flow, sigma, r_max): - """ - Compute Gaussian density map from optical flow using CUDA acceleration. - - This function uses a custom CUDA kernel for parallel computation of density maps, - providing significant speedup over CPU implementation for large flow fields. - - Args: - flow: Optical flow field with shape (height, width, 2), can be numpy array - or OpenCV CUDA GpuMat - sigma: Gaussian kernel standard deviation - r_max: Maximum radius for Gaussian contribution (cutoff distance) - - Returns: - np.ndarray: Density map with shape (height, width), transferred from GPU to CPU - - Note: - Requires CuPy and CUDA-capable GPU. Falls back to CPU if unavailable. - """ - # Precompute 2*sigma^2 for efficiency - twoSigma2 = cp.float32(2.0 * sigma**2) - r_max = cp.float32(r_max) - - # Convert input to CuPy arrays - # Handle both numpy arrays and OpenCV CUDA GpuMat - if type(flow).__module__ == np.__name__: - # NumPy array: convert directly to CuPy - flowx = cp.array(flow[:, :, 0], dtype=cp.float32) - flowy = cp.array(flow[:, :, 1], dtype=cp.float32) - else: - # OpenCV CUDA GpuMat: use zero-copy interface - flowx, flowy = cv2.cuda.split(flow) - flowx = cp.array(CudaArrayInterface(flowx), dtype=cp.float32, copy=True) - flowy = cp.array(CudaArrayInterface(flowy), dtype=cp.float32, copy=True) - - # Initialize output array on GPU - result = cp.zeros((flowx.shape[0], flowx.shape[1]), dtype=cp.float32) - - # Configure CUDA kernel launch parameters - block_dim = (flowx.shape[1], ) # Threads per block (one thread per column) - grid = (flowx.shape[0], ) # Number of blocks (one block per row) - - # Launch CUDA kernel - funcupy(grid, block_dim, (flowx, flowy, result, - cp.int32(flowx.shape[1]), - cp.int32(flowx.shape[0]), - twoSigma2, r_max)) - - # Transfer result from GPU to CPU - return result.get() diff --git a/src/orisys/nHHD.py b/src/orisys/nHHD.py deleted file mode 100644 index 85d3c4f..0000000 --- a/src/orisys/nHHD.py +++ /dev/null @@ -1,798 +0,0 @@ -''' -Copyright (c) 2015, Harsh Bhatia (bhatia4@llnl.gov) -All rights reserved. - -Redistribution and use in source and binary forms, with or without -modification, are permitted provided that the following conditions are met: -1. Redistributions of source code must retain the above copyright notice, this -list of conditions and the following disclaimer. -2. Redistributions in binary form must reproduce the above copyright notice, -this list of conditions and the following disclaimer in the documentation -and/or other materials provided with the distribution. -THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND -ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED -WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE -DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE -FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL -DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER -CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, -OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -''' - -import numpy -from .nPoisson import PoissonSolver -from .StructuredGrid import StructuredGrid -from .UnstructuredGrid import UnstructuredGrid -# import matplotlib.pyplot as plt - -# Try to import CuPy for GPU acceleration -try: - import cupy as cp - CUPY_AVAILABLE = True - - # Try to import CuPy streams for parallel execution - try: - from cupy.cuda import Stream - STREAMS_AVAILABLE = True - except ImportError: - Stream = None - STREAMS_AVAILABLE = False - - # Create fused elementwise kernel for torque computation - # Computes cross products and combines them in a single kernel pass - try: - _torque_kernel = cp.ElementwiseKernel( - 'float32 r_x, float32 r_y, float32 grid_0, float32 grid_1, ' - 'float32 max_0, float32 max_1, float32 min_0, float32 min_1', - 'float32 torque', - ''' - float cross_p = (grid_0 - max_0) * r_y - (grid_1 - max_1) * r_x; - float cross_n = (grid_0 - min_0) * r_y - (grid_1 - min_1) * r_x; - torque = cross_p + cross_n; - ''', - 'torque_kernel' - ) - KERNEL_AVAILABLE = True - except Exception: - _torque_kernel = None - KERNEL_AVAILABLE = False - -except ImportError: - cp = None - CUPY_AVAILABLE = False - STREAMS_AVAILABLE = False - KERNEL_AVAILABLE = False - Stream = None - _torque_kernel = None - -# --------------------------------------------------------------------------- -class nHHD(object): - - # constructor - def __init__(self, **kwargs): - - args = kwargs.keys() - - if ('grid' in args) == ('points' in args): - raise SyntaxError("nHHD Solver needs either shape of a regular grid or the points in an rectilinear/unstructured grid") - - if ('grid' in args) != ('spacings' in args): - raise SyntaxError("nHHD needs spacings for regular grid") - - if 'grid' in args: - - self.mesh = StructuredGrid(grid=kwargs['grid'], spacings=kwargs['spacings']) - self.psolver = PoissonSolver(grid=self.mesh.dims, spacings=self.mesh.dx) - - elif 'points' in args: - - if 'simplices' in args: - self.mesh = UnstructuredGrid(vertices=kwargs['points'], simplices=kwargs['simplices']) - else: - self.mesh = UnstructuredGrid(vertices=kwargs['points']) - - self.psolver = PoissonSolver(points=self.mesh.vertices, pvolumes=self.mesh.pvolumes) - - self.dim = self.mesh.dim - if (self.dim != 2) and (self.dim != 3): - raise ValueError("nHHD Solver works only for 2 and 3 dimensions") - - self.psolver.prepare(False) - - # create the 3-component decomposition - def decompose(self, vfield, verbose=0): - - if vfield.shape[-1] != self.dim: - raise ValueError("nHHD.decompose requires a valid-dimensional vector field") - - if verbose > 0: - print( '\n ==> nhhd.decompose', vfield.shape) - - if self.dim == 2: - - # compute div curl - (self.div, self.curlw) = self.mesh.divcurl(vfield) - - if verbose > 0: - print( 'div =', self.div.shape, self.div.min(), self.div.max()) - print( 'curlw =', self.curlw.shape, self.curlw.min(), self.curlw.max()) - - # compute potentials - self.nD = self.psolver.solve_gpu(self.div, verbose > 0) #1) - if verbose > 0: - print( 'nD =', self.nD.shape, self.nD.min(), self.nD.max()) - - self.nRu = self.psolver.solve_gpu(self.curlw, verbose > 0) #1) - if verbose > 0: - print( 'nR =', self.nRu.shape, self.nRu.min(), self.nRu.max()) - - # compute fields as gradients of potentials - self.d = self.mesh.gradient(self.nD) - if verbose > 0: - print( 'd =', self.d.shape, self.d.min(), self.d.max()) - - self.r = self.mesh.rotated_gradient(self.nRu) - if verbose > 0: - print( 'r =', self.r.shape, self.r.min(), self.r.max()) - else: - - # compute div curl - (self.div, self.curlu, self.curlv, self.curlw) = self.mesh.divcurl(vfield) - - if verbose > 0: - print( 'div =', self.div.shape, self.div.min(), self.div.max()) - print( 'curlu =', self.curlu.shape, self.curlu.min(), self.curlu.max()) - print( 'curlv =', self.curlv.shape, self.curlv.min(), self.curlv.max()) - print( 'curlw =', self.curlw.shape, self.curlw.min(), self.curlw.max()) - - # compute potentials - self.nD = self.psolver.solve_gpu(self.div, verbose > 0) #1) - if verbose > 0: - print( 'nD =', self.nD.shape, self.nD.min(), self.nD.max()) - - self.nRu = self.psolver.solve_gpu(self.curlu, verbose > 0) - self.nRu *= -1.0 - if verbose > 0: - print( 'nRu =', self.nRu.shape, self.nRu.min(), self.nRu.max()) - - self.nRv = self.psolver.solve_gpu(self.curlv, verbose > 0) - self.nRv *= -1.0 - if verbose > 0: - print( 'nRv =', self.nRv.shape, self.nRv.min(), self.nRv.max()) - - self.nRw = self.psolver.solve_gpu(self.curlw, verbose > 0) - self.nRw *= -1.0 - - if verbose > 0: - print( 'nRw =', self.nRw.shape, self.nRw.min(), self.nRw.max()) - - # compute divergent field as gradient of D potential - self.d = self.mesh.gradient(self.nD) - if verbose > 0: - print( 'd =', self.d.shape, self.d.min(), self.d.max()) - - # temporary use of r - self.r = numpy.stack((self.nRu, self.nRv, self.nRw), axis = -1) - - # compute rotational field as curl of R potential - (cu, cv, cw) = self.mesh.curl3D(self.r) - self.r = numpy.stack((cu, cv, cw), axis = -1) - if verbose > 0: - print( 'r =', self.r.shape, self.r.min(), self.r.max()) - - self.h = numpy.add(self.d, self.r) - numpy.subtract(vfield, self.h, self.h) - - if verbose > 0: - print( 'h =', self.h.shape, self.h.min(), self.h.max()) - - def decompose_gpu(self, vfield, verbose=0): - """ - GPU-optimized decomposition that keeps intermediate results on GPU. - Returns GPU arrays to avoid CPU-GPU transfers. - - Args: - vfield: Vector field (numpy array or CuPy array) - verbose: Verbosity level (default: 0) - - Returns: - Tuple of GPU arrays: (d_gpu, h_gpu, r_gpu, nRu_gpu) for 2D fields - """ - if not CUPY_AVAILABLE: - if verbose > 0: - print("CuPy not available, falling back to CPU decompose") - self.decompose(vfield, verbose) - return None - - if vfield.shape[-1] != self.dim: - raise ValueError("nHHD.decompose requires a valid-dimensional vector field") - - if verbose > 0: - print('\n ==> nhhd.decompose_gpu', vfield.shape) - - # Convert input to GPU if needed - if isinstance(vfield, numpy.ndarray): - vfield_gpu = cp.asarray(vfield, dtype=cp.float32) - else: - vfield_gpu = vfield - - if self.dim == 2: - # Compute div curl on GPU - dx0, dx1 = self.mesh.dx[0], self.mesh.dx[1] - - dudy, dudx = cp.gradient(vfield_gpu[:,:,0], dx0, dx1) - dvdy, dvdx = cp.gradient(vfield_gpu[:,:,1], dx0, dx1) - - div_gpu = dudx + dvdy - curlw_gpu = dvdx - dudy - - if verbose > 0: - print('div =', div_gpu.shape, float(div_gpu.min()), float(div_gpu.max())) - print('curlw =', curlw_gpu.shape, float(curlw_gpu.min()), float(curlw_gpu.max())) - - # Compute potentials on GPU (solve_gpu returns numpy, so we convert back) - nD_cpu = self.psolver.solve_gpu(cp.asnumpy(div_gpu), verbose > 0) - nRu_cpu = self.psolver.solve_gpu(cp.asnumpy(curlw_gpu), verbose > 0) - - # Keep potentials on GPU - nD_gpu = cp.asarray(nD_cpu, dtype=cp.float32) - nRu_gpu = cp.asarray(nRu_cpu, dtype=cp.float32) - - if verbose > 0: - print('nD =', nD_gpu.shape, float(nD_gpu.min()), float(nD_gpu.max())) - print('nR =', nRu_gpu.shape, float(nRu_gpu.min()), float(nRu_gpu.max())) - - # Compute gradients on GPU - ddy, ddx = cp.gradient(nD_gpu, dx0, dx1) - d_gpu = cp.stack((ddx, ddy), axis=-1) - - if verbose > 0: - print('d =', d_gpu.shape, float(d_gpu.min()), float(d_gpu.max())) - - # Compute rotated gradient on GPU - rdy, rdx = cp.gradient(nRu_gpu, dx0, dx1) - rdy *= -1.0 - r_gpu = cp.stack((rdy, rdx), axis=-1) - - if verbose > 0: - print('r =', r_gpu.shape, float(r_gpu.min()), float(r_gpu.max())) - - # Compute h on GPU - h_gpu = d_gpu + r_gpu - h_gpu = vfield_gpu - h_gpu - - if verbose > 0: - print('h =', h_gpu.shape, float(h_gpu.min()), float(h_gpu.max())) - - # Store GPU arrays - self._d_gpu = d_gpu - self._h_gpu = h_gpu - self._r_gpu = r_gpu - self._nRu_gpu = nRu_gpu - - # Also store CPU copies for compatibility - self.d = cp.asnumpy(d_gpu) - self.h = cp.asnumpy(h_gpu) - self.r = cp.asnumpy(r_gpu) - self.nRu = cp.asnumpy(nRu_gpu) - self.nD = nD_cpu - - return (d_gpu, h_gpu, r_gpu, nRu_gpu) - - else: - # For 3D, fall back to CPU for now - if verbose > 0: - print("3D GPU decompose not yet implemented, using CPU") - self.decompose(vfield, verbose) - return None - - def decompose_batch(self, vfield, verbose=0): - - if vfield.shape[-1] != self.dim: - raise ValueError("nHHD.decompose requires a valid-dimensional vector field") - - if verbose > 0: - print( '\n ==> nhhd.decompose', vfield.shape) - - if self.dim == 2: - - # compute div curl using GPU - (self.div, self.curlw) = self.mesh.divcurl(vfield) - - if verbose > 0: - print( 'div =', self.div.shape, self.div.min(), self.div.max()) - print( 'curlw =', self.curlw.shape, self.curlw.min(), self.curlw.max()) - - # compute potentials using batch GPU solver - batch_inputs = [self.div, self.curlw] - batch_results = self.psolver.solve_gpu_batch(batch_inputs, verbose > 0) - - self.nD = batch_results[0] - self.nRu = batch_results[1] - - if verbose > 0: - print( 'nD =', self.nD.shape, self.nD.min(), self.nD.max()) - print( 'nR =', self.nRu.shape, self.nRu.min(), self.nRu.max()) - - # compute fields as gradients of potentials - self.d = self.mesh.gradient(self.nD) - if verbose > 0: - print( 'd =', self.d.shape, self.d.min(), self.d.max()) - - self.r = self.mesh.rotated_gradient(self.nRu) - if verbose > 0: - print( 'r =', self.r.shape, self.r.min(), self.r.max()) - else: - - # compute div curl - (self.div, self.curlu, self.curlv, self.curlw) = self.mesh.divcurl(vfield) - - if verbose > 0: - print( 'div =', self.div.shape, self.div.min(), self.div.max()) - print( 'curlu =', self.curlu.shape, self.curlu.min(), self.curlu.max()) - print( 'curlv =', self.curlv.shape, self.curlv.min(), self.curlv.max()) - print( 'curlw =', self.curlw.shape, self.curlw.min(), self.curlw.max()) - - # compute potentials using batch GPU solver - batch_inputs = [self.div, self.curlu, self.curlv, self.curlw] - batch_results = self.psolver.solve_gpu_batch(batch_inputs, verbose > 0) - - self.nD = batch_results[0] - self.nRu = batch_results[1] * -1.0 - self.nRv = batch_results[2] * -1.0 - self.nRw = batch_results[3] * -1.0 - - if verbose > 0: - print( 'nD =', self.nD.shape, self.nD.min(), self.nD.max()) - print( 'nRu =', self.nRu.shape, self.nRu.min(), self.nRu.max()) - print( 'nRv =', self.nRv.shape, self.nRv.min(), self.nRv.max()) - print( 'nRw =', self.nRw.shape, self.nRw.min(), self.nRw.max()) - - # compute divergent field as gradient of D potential - self.d = self.mesh.gradient(self.nD) - if verbose > 0: - print( 'd =', self.d.shape, self.d.min(), self.d.max()) - - # temporary use of r - self.r = numpy.stack((self.nRu, self.nRv, self.nRw), axis = -1) - - # compute rotational field as curl of R potential - (cu, cv, cw) = self.mesh.curl3D(self.r) - self.r = numpy.stack((cu, cv, cw), axis = -1) - if verbose > 0: - print( 'r =', self.r.shape, self.r.min(), self.r.max()) - - self.h = numpy.add(self.d, self.r) - numpy.subtract(vfield, self.h, self.h) - - if verbose > 0: - print( 'h =', self.h.shape, self.h.min(), self.h.max()) - - def get_force(self, vfield, verbose=False): - - mgrid = numpy.mgrid[0: vfield.shape[0], 0: vfield.shape[1]] - grid = numpy.stack((mgrid[0], mgrid[1]), axis =2) - # print(grid.shape) - - - self.decompose(vfield, verbose) - - self.Sn = numpy.sum(numpy.absolute(self.d)) - # print("Sn is:", self.Sn) - - # self.Stx = abs(numpy.sum(self.h[:, :, 0])) - # self.Sty = abs(numpy.sum(self.h[:, :, 1])) - self.Stx = numpy.sum(self.h[:, :, 0]) - self.Sty = numpy.sum(self.h[:, :, 1]) - # print("Stx is:", self.Stx) - # print("Sty is:", self.Sty) - - nR_max_loc = numpy.unravel_index(numpy.argmax(self.nRu, axis = None), self.nRu.shape) - nR_min_loc = numpy.unravel_index(numpy.argmin(self.nRu, axis = None), self.nRu.shape) - - # print(nR_max_loc) - - r_p = grid - nR_max_loc - r_n = grid - nR_min_loc - - self.Stor = numpy.sum(numpy.cross(r_p, self.r)) + numpy.sum(numpy.cross(r_n, self.r)) - # print("Stor is:", self.Stor) - - output = numpy.array([self.Sn, self.Stx, self.Sty, self.Stor]) - - flag = False - if self.Sn >=300 and self.Stx >= 300 and self.Sty >= 300 and numpy.abs(self.Stor) >= 1000: - flag = True - # print(flag) - return flag,output - # return self.d, self.h, numpy.cross(r_p, self.r) , numpy.cross(r_n, self.r) - - def _init_force_gpu_cache(self, shape): - """ - Initialize cached GPU arrays for force computation. - Call once when shape is known to avoid repeated allocations. - - Args: - shape: Tuple (height, width) of the vector field - """ - if not CUPY_AVAILABLE: - return - - # Cache grid arrays (only need to create once per shape) - self._force_gpu_cache = { - 'shape': shape, - 'mgrid_0': cp.arange(shape[0], dtype=cp.float32)[:, None], - 'mgrid_1': cp.arange(shape[1], dtype=cp.float32)[None, :], - } - # Pre-broadcast for efficiency - self._force_gpu_cache['mgrid_0_bc'] = cp.broadcast_to( - self._force_gpu_cache['mgrid_0'], shape - ).copy() # .copy() to make contiguous - self._force_gpu_cache['mgrid_1_bc'] = cp.broadcast_to( - self._force_gpu_cache['mgrid_1'], shape - ).copy() - - def get_force_gpu(self, vfield, verbose=False): - """ - GPU-accelerated force computation using CuPy. - Falls back to CPU if CuPy is not available. - - Args: - vfield: Vector field (numpy array) - verbose: Print debug information (default: False) - - Returns: - tuple: (flag, output) where: - - flag: Boolean indicating if force thresholds are met - - output: Array [Sn, Stx, Sty, Stor] containing force components - """ - if not CUPY_AVAILABLE: - if verbose: - print("CuPy not available, falling back to CPU get_force") - return self.get_force(vfield, verbose) - - # Decompose the vector field (uses GPU internally) - self.decompose(vfield, verbose) - - # Initialize cache if needed or shape changed - shape = (vfield.shape[0], vfield.shape[1]) - if not hasattr(self, '_force_gpu_cache') or self._force_gpu_cache is None or self._force_gpu_cache['shape'] != shape: - self._init_force_gpu_cache(shape) - - # Transfer data to GPU (consider keeping on GPU in decompose for further optimization) - d_gpu = cp.asarray(self.d, dtype=cp.float32) - h_gpu = cp.asarray(self.h, dtype=cp.float32) - r_gpu = cp.asarray(self.r, dtype=cp.float32) - nRu_gpu = cp.asarray(self.nRu, dtype=cp.float32) - - # Compute all sums in parallel - fused for better performance - # Sn: sum of absolute values of divergent field (both components) - Sn_gpu = cp.sum(cp.abs(d_gpu)) - - # Stx and Sty: sum of harmonic field components - Stx_gpu = cp.sum(h_gpu[:, :, 0]) - Sty_gpu = cp.sum(h_gpu[:, :, 1]) - - # Find max/min locations of nRu using single pass - nR_max_idx = int(cp.argmax(nRu_gpu)) - nR_min_idx = int(cp.argmin(nRu_gpu)) - - # Convert flat indices to 2D coordinates (on CPU - faster for single values) - nR_max_loc = numpy.unravel_index(nR_max_idx, nRu_gpu.shape) - nR_min_loc = numpy.unravel_index(nR_min_idx, nRu_gpu.shape) - - # Use cached grid arrays - mgrid_0 = self._force_gpu_cache['mgrid_0_bc'] - mgrid_1 = self._force_gpu_cache['mgrid_1_bc'] - - # Compute torque using fused operations - # r_p and r_n displacements + cross products in one pass - max_0, max_1 = float(nR_max_loc[0]), float(nR_max_loc[1]) - min_0, min_1 = float(nR_min_loc[0]), float(nR_min_loc[1]) - - r_x = r_gpu[:, :, 0] - r_y = r_gpu[:, :, 1] - - # Fused cross product computation - # cross_p = (mgrid_0 - max_0) * r_y - (mgrid_1 - max_1) * r_x - # cross_n = (mgrid_0 - min_0) * r_y - (mgrid_1 - min_1) * r_x - Stor_gpu = cp.sum( - (mgrid_0 - max_0) * r_y - (mgrid_1 - max_1) * r_x + - (mgrid_0 - min_0) * r_y - (mgrid_1 - min_1) * r_x - ) - - # Single batch transfer to CPU (more efficient than 4 separate transfers) - results = cp.asnumpy(cp.array([Sn_gpu, Stx_gpu, Sty_gpu, Stor_gpu])) - - self.Sn = float(results[0]) - self.Stx = float(results[1]) - self.Sty = float(results[2]) - self.Stor = float(results[3]) - - output = results - - # Check thresholds - flag = (self.Sn >= 300 and self.Stx >= 300 and self.Sty >= 300 and numpy.abs(self.Stor) >= 1000) - - return flag, output - - def get_force_gpu_v2(self, vfield, verbose=False): - """ - Further optimized GPU force computation using CuPy ElementwiseKernel. - Fuses all operations into a single kernel for maximum performance. - - Args: - vfield: Vector field (numpy array) - verbose: Print debug information (default: False) - - Returns: - tuple: (flag, output) where: - - flag: Boolean indicating if force thresholds are met - - output: Array [Sn, Stx, Sty, Stor] containing force components - """ - if not CUPY_AVAILABLE: - if verbose: - print("CuPy not available, falling back to CPU get_force") - return self.get_force(vfield, verbose) - - # Decompose the vector field - self.decompose(vfield, verbose) - - # Keep data on GPU - avoid transfers by computing directly - d_gpu = cp.asarray(self.d, dtype=cp.float32) - h_gpu = cp.asarray(self.h, dtype=cp.float32) - r_gpu = cp.asarray(self.r, dtype=cp.float32) - nRu_gpu = cp.asarray(self.nRu, dtype=cp.float32) - - # Use reduction kernels for better performance - # Compute Sn using L1 norm (sum of absolute values) - Sn = float(cp.linalg.norm(d_gpu.ravel(), ord=1)) - - # Compute Stx, Sty - Stx = float(cp.sum(h_gpu[:, :, 0])) - Sty = float(cp.sum(h_gpu[:, :, 1])) - - # Find extrema indices - nR_max_idx = int(cp.argmax(nRu_gpu)) - nR_min_idx = int(cp.argmin(nRu_gpu)) - - shape = nRu_gpu.shape - max_0, max_1 = nR_max_idx // shape[1], nR_max_idx % shape[1] - min_0, min_1 = nR_min_idx // shape[1], nR_min_idx % shape[1] - - # Initialize cache if needed - if not hasattr(self, '_force_gpu_cache') or self._force_gpu_cache is None or self._force_gpu_cache['shape'] != shape: - self._init_force_gpu_cache(shape) - - mgrid_0 = self._force_gpu_cache['mgrid_0_bc'] - mgrid_1 = self._force_gpu_cache['mgrid_1_bc'] - - r_x = r_gpu[:, :, 0] - r_y = r_gpu[:, :, 1] - - # Compute torque with fused sum - Stor = float(cp.sum( - (mgrid_0 - max_0) * r_y - (mgrid_1 - max_1) * r_x + - (mgrid_0 - min_0) * r_y - (mgrid_1 - min_1) * r_x - )) - - self.Sn, self.Stx, self.Sty, self.Stor = Sn, Stx, Sty, Stor - output = numpy.array([Sn, Stx, Sty, Stor], dtype=numpy.float32) - - flag = (Sn >= 300 and Stx >= 300 and Sty >= 300 and abs(Stor) >= 1000) - - return flag, output - - def clear_force_gpu_cache(self): - """Clear the GPU cache for force computation to free memory.""" - if hasattr(self, '_force_gpu_cache'): - self._force_gpu_cache = None - if CUPY_AVAILABLE: - cp.get_default_memory_pool().free_all_blocks() - - def get_force_gpu_stream(self, vfield, verbose=False): - """ - Fully optimized GPU force computation using: - - decompose_gpu to keep data on GPU throughout - - CUDA streams for parallel reductions - - Fused kernel operations - - Args: - vfield: Vector field (numpy array) - verbose: Print debug information (default: False) - - Returns: - tuple: (flag, output) where: - - flag: Boolean indicating if force thresholds are met - - output: Array [Sn, Stx, Sty, Stor] containing force components - """ - if not CUPY_AVAILABLE: - if verbose: - print("CuPy not available, falling back to CPU get_force") - return self.get_force(vfield, verbose) - - # Use GPU decomposition to keep data on GPU - gpu_result = self.decompose_gpu(vfield, verbose) - - if gpu_result is None: - # Fallback for 3D or if decompose_gpu failed - return self.get_force_gpu(vfield, verbose) - - d_gpu, h_gpu, r_gpu, nRu_gpu = gpu_result - - # Initialize cache if needed or shape changed - shape = (vfield.shape[0], vfield.shape[1]) - if not hasattr(self, '_force_gpu_cache') or self._force_gpu_cache is None or self._force_gpu_cache['shape'] != shape: - self._init_force_gpu_cache(shape) - - # Use CUDA streams for parallel computation if available - if STREAMS_AVAILABLE: - stream1 = Stream(non_blocking=True) - stream2 = Stream(non_blocking=True) - stream3 = Stream(non_blocking=True) - - # Stream 1: Compute Sn - with stream1: - Sn_gpu = cp.sum(cp.abs(d_gpu)) - - # Stream 2: Compute Stx, Sty - with stream2: - Stx_gpu = cp.sum(h_gpu[:, :, 0]) - Sty_gpu = cp.sum(h_gpu[:, :, 1]) - - # Stream 3: Find extrema and compute torque - with stream3: - nR_max_idx = int(cp.argmax(nRu_gpu)) - nR_min_idx = int(cp.argmin(nRu_gpu)) - - max_0, max_1 = nR_max_idx // shape[1], nR_max_idx % shape[1] - min_0, min_1 = nR_min_idx // shape[1], nR_min_idx % shape[1] - - mgrid_0 = self._force_gpu_cache['mgrid_0_bc'] - mgrid_1 = self._force_gpu_cache['mgrid_1_bc'] - - r_x = r_gpu[:, :, 0] - r_y = r_gpu[:, :, 1] - - Stor_gpu = cp.sum( - (mgrid_0 - max_0) * r_y - (mgrid_1 - max_1) * r_x + - (mgrid_0 - min_0) * r_y - (mgrid_1 - min_1) * r_x - ) - - # Synchronize all streams - stream1.synchronize() - stream2.synchronize() - stream3.synchronize() - else: - # Without streams, compute sequentially but still on GPU - Sn_gpu = cp.sum(cp.abs(d_gpu)) - Stx_gpu = cp.sum(h_gpu[:, :, 0]) - Sty_gpu = cp.sum(h_gpu[:, :, 1]) - - nR_max_idx = int(cp.argmax(nRu_gpu)) - nR_min_idx = int(cp.argmin(nRu_gpu)) - - max_0, max_1 = nR_max_idx // shape[1], nR_max_idx % shape[1] - min_0, min_1 = nR_min_idx // shape[1], nR_min_idx % shape[1] - - mgrid_0 = self._force_gpu_cache['mgrid_0_bc'] - mgrid_1 = self._force_gpu_cache['mgrid_1_bc'] - - r_x = r_gpu[:, :, 0] - r_y = r_gpu[:, :, 1] - - Stor_gpu = cp.sum( - (mgrid_0 - max_0) * r_y - (mgrid_1 - max_1) * r_x + - (mgrid_0 - min_0) * r_y - (mgrid_1 - min_1) * r_x - ) - - # Single batch transfer to CPU - results = cp.asnumpy(cp.array([Sn_gpu, Stx_gpu, Sty_gpu, Stor_gpu])) - - self.Sn = float(results[0]) - self.Stx = float(results[1]) - self.Sty = float(results[2]) - self.Stor = float(results[3]) - - output = results - - # Check thresholds - flag = (self.Sn >= 300 and self.Stx >= 300 and self.Sty >= 300 and numpy.abs(self.Stor) >= 1000) - - return flag, output - - def get_force_gpu_kernel(self, vfield, verbose=False): - """ - Maximum performance GPU force computation using: - - decompose_gpu to keep data on GPU throughout - - Custom CUDA kernel for fused torque computation - - Minimized memory allocations - - Args: - vfield: Vector field (numpy array) - verbose: Print debug information (default: False) - - Returns: - tuple: (flag, output) where: - - flag: Boolean indicating if force thresholds are met - - output: Array [Sn, Stx, Sty, Stor] containing force components - """ - if not CUPY_AVAILABLE: - if verbose: - print("CuPy not available, falling back to CPU get_force") - return self.get_force(vfield, verbose) - - # Use GPU decomposition to keep data on GPU - gpu_result = self.decompose_gpu(vfield, verbose) - - if gpu_result is None: - return self.get_force_gpu(vfield, verbose) - - d_gpu, h_gpu, r_gpu, nRu_gpu = gpu_result - - # Initialize cache if needed - shape = (vfield.shape[0], vfield.shape[1]) - if not hasattr(self, '_force_gpu_cache') or self._force_gpu_cache is None or self._force_gpu_cache['shape'] != shape: - self._init_force_gpu_cache(shape) - - # Compute Sn using L1 norm (faster than sum(abs())) - Sn = float(cp.linalg.norm(d_gpu.ravel(), ord=1)) - - # Compute Stx, Sty - Stx = float(cp.sum(h_gpu[:, :, 0])) - Sty = float(cp.sum(h_gpu[:, :, 1])) - - # Find extrema - nR_max_idx = int(cp.argmax(nRu_gpu)) - nR_min_idx = int(cp.argmin(nRu_gpu)) - - max_0 = cp.float32(nR_max_idx // shape[1]) - max_1 = cp.float32(nR_max_idx % shape[1]) - min_0 = cp.float32(nR_min_idx // shape[1]) - min_1 = cp.float32(nR_min_idx % shape[1]) - - mgrid_0 = self._force_gpu_cache['mgrid_0_bc'] - mgrid_1 = self._force_gpu_cache['mgrid_1_bc'] - - r_x = r_gpu[:, :, 0] - r_y = r_gpu[:, :, 1] - - # Use custom kernel if available for fused torque computation - if KERNEL_AVAILABLE and _torque_kernel is not None: - torque_vals = _torque_kernel(r_x, r_y, mgrid_0, mgrid_1, max_0, max_1, min_0, min_1) - Stor = float(cp.sum(torque_vals)) - else: - # Fallback to standard computation - Stor = float(cp.sum( - (mgrid_0 - max_0) * r_y - (mgrid_1 - max_1) * r_x + - (mgrid_0 - min_0) * r_y - (mgrid_1 - min_1) * r_x - )) - - self.Sn, self.Stx, self.Sty, self.Stor = Sn, Stx, Sty, Stor - output = numpy.array([Sn, Stx, Sty, Stor], dtype=numpy.float32) - - flag = (Sn >= 300 and Stx >= 300 and Sty >= 300 and abs(Stor) >= 1000) - - return flag, output - - def get_shear(self): - - return self.h - - def get_normal(self): - - return self.d - - # def plot_2d_information(self,value): - # fig,ax = plt.subplots() - # x,y = [],[] - # for i in range (1000000000000): - # x.append(i) - # y.append(value) - # ax.cla() - # ax.plot(x,y,'b',lw=1) - - # plt.pause(0.1) -# --------------------------------------------------------------------------- diff --git a/src/orisys/nPoisson.py b/src/orisys/nPoisson.py deleted file mode 100644 index 4bcee3b..0000000 --- a/src/orisys/nPoisson.py +++ /dev/null @@ -1,518 +0,0 @@ -''' -Copyright (c) 2015, Harsh Bhatia (bhatia4@llnl.gov) -All rights reserved. - -Redistribution and use in source and binary forms, with or without -modification, are permitted provided that the following conditions are met: -1. Redistributions of source code must retain the above copyright notice, this -list of conditions and the following disclaimer. -2. Redistributions in binary form must reproduce the above copyright notice, -this list of conditions and the following disclaimer in the documentation -and/or other materials provided with the distribution. -THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND -ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED -WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE -DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE -FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL -DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER -CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, -OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -''' - -import numpy -from scipy import signal, spatial -from .timer import Timer - -# Try to import CuPy for GPU acceleration -try: - import cupy as cp - import cupyx.scipy.signal as cp_signal - CUPY_AVAILABLE = True -except ImportError: - cp = None - cp_signal = None - CUPY_AVAILABLE = False - - -''' -# ------------------------------------------------------------------------------ - Poisson Solver for 1D, 2D, and 3D scalar fields. - Data may be defined on a rectlinear grid or unstructured points. - - 1) create an object - PoissonSolver(args) - 'solver' : type of solver - 'S' for spatial convolution - only option for unstructured points - 'F' for convolution in frequency space - default option for regular grids - - --> for regular grid - 'grid' : shape of the rectlinear grid - tuple of dimensions (Z,Y,X) - 'spacings' : spacings of the rectlinear grid - tuple of spacings (dz,dy,dx) - - --> for unstrucutred points - 'points' : unstructured points - numpy array of shape (N, dim) - 'pvolumes' : point volumes for each point - numpy array of shape (N, 1) - - 2) call prepare() - this initializes internal datastructures depending upon the mesh. - - Once initialized, solve() can be called any number of times to - solve the Poisson equation on the same domain. If the domain - changes (e.g., a Lagrangian mesh), create a new object. - - 3) call solve(f), where, f is a numpy array defining a scalar field - - for rectlinear grid, f.shape = grid.shape - for unstructured points, f.shape = (points.shape[0], 1) - -# ------------------------------------------------------------------------------ -''' - -class PoissonSolver(object): - - # constructor - def __init__(self, **kwargs): - - args = kwargs.keys() - - if ('grid' in args) == ('points' in args): - raise SyntaxError("Poisson Solver needs either shape of a regular grid or the points in an rectilinear/unstructured grid") - - if ('grid' in args) != ('spacings' in args): - raise SyntaxError("Poisson Solver needs spacings for regular grid") - - if ('points' in args) != ('pvolumes' in args): - raise SyntaxError("Poisson Solver needs point volumes for unstructured points") - - # (default) data type - self.dtype = numpy.float32 - if 'dtype' in args: - self.dtype = kwargs['dtype'] - - if 'grid' in args: - self.ptype = 'G' - self.stype = 'F' # default - self.gdims = kwargs['grid'] - self.gdx = kwargs['spacings'] - self.dim = len(self.gdims) - if self.dim != len(self.gdx): - raise ValueError("Dimensions of spacings should match that of the grid") - - elif 'points' in args: - self.ptype = 'P' - self.stype = 'S' # only option - self.points = kwargs['points'] - self.pvolumes = kwargs['pvolumes'] - self.dim = len(self.points.shape) - if self.pvolumes.shape[0] != self.points.shape[0]: - raise ValueError("Number of pvolumes should match that of the points") - - if 'solver' in args: - if (kwargs['solver'] != 'S' and kwargs['solver'] != 'F'): - raise SyntaxError('Solver can be only of type S(patial) or F(requency)') - - if (kwargs['solver'] == 'F') and ('points' in args): - raise SyntaxError('Only Spatial solver supported for unstructured points') - - if (kwargs['solver'] == 'S') and ('grid' in args): - raise SyntaxError('Frequncy solver strongly suggested for structured grids') - - self.stype = kwargs['solver'] - - if (self.dim > 3): - raise ValueError("Poisson Solver works only for 1, 2, or 3 dimensions") - - self.ready = False - self._gpu_cache = None # Cache for GPU arrays - # if (self.ptype == 'G'): - # print( 'PoissonSolver:', self.dim, 'D grid =', self.gdims, 'with spacings', self.gdx) - # elif (self.ptype == 'P'): - # print( 'PoissonSolver', self.dim, 'D points =', self.points.shape) - - # -------------------------------------------------------------------------- - # create the Green's function in appropriate representation - # -------------------------------------------------------------------------- - - # generate a grid -(N-1) to (N-1) to define distance function - # where N is the size of the original grid - def _create_radialGrid(self, verbose): - - mtimer = Timer() - if verbose: - print( ' - creating distance kernel:', end ='') - - rdims = tuple(2*d - 1 for d in self.gdims) - - R = numpy.indices(rdims, dtype=self.dtype) - for d in range(self.dim): - R[d] += 1-self.gdims[d] - R[d] *= self.gdx[d] - - if self.dim == 1: - numpy.absolute(R[0], R[0]) - R = R[0] - - elif self.dim == 2: - numpy.hypot(R[1], R[0], R[0]) - R = R[0] - - # half of the average grid spacing - zval = 0.5 * (self.gdx[0]+self.gdx[1]) / 2.0 - R[self.gdims[0]-1, self.gdims[1]-1] = zval - - elif self.dim == 3: - numpy.hypot(R[2], R[1], R[1]) - numpy.hypot(R[1], R[0], R[0]) - R = R[0] - - # half of the average grid spacing - zval = 0.5 * (self.gdx[0]+self.gdx[1]+self.gdx[2]) / 3.0 - R[self.gdims[0]-1, self.gdims[1]-1, self.gdims[2]-1] = zval - - if verbose: - print( R.shape, R.min(), R.max()) - mtimer.end() - - return R - - # -------------------------------------------------------------------------- - # scale a function by point volumes for correct integration - def _scale(self, func): - - if self.ptype == 'P': - return numpy.multiply(func, self.pvolumes) - - elif self.ptype == 'G': - sfunc = func * numpy.prod(self.gdx) - - ''' - # volume of the boundary elements is only half the voxel - if self.dim == 1: - sfunc[ 0] *= 0.5 - sfunc[-1] *= 0.5 - elif self.dim == 2: - sfunc[ 0, :] *= 0.5 - sfunc[-1, :] *= 0.5 - sfunc[ :, 0] *= 0.5 - sfunc[ :,-1] *= 0.5 - elif self.dim == 3: - sfunc[ 0, :, :] *= 0.5 - sfunc[-1, :, :] *= 0.5 - sfunc[ :, 0, :] *= 0.5 - sfunc[ :,-1, :] *= 0.5 - sfunc[ :, :, 0] *= 0.5 - sfunc[ :, :,-1] *= 0.5 - ''' - return sfunc - - def prepare(self, verbose=False): - - # ---------------------------------------------------------------------- - gtimer = Timer() - if verbose: - print( '\nInitializing Poisson solver, type =', self.stype) - - # ---------------------------------------------------------------------- - # generate a grid to store pairwise distances to store Green's function - self.G = None # definition depends upon solver type - - # ---------------------------------------------------------------------- - # for Fourier solution, - if self.stype == 'F': - - if self.ptype == 'G': - self.G = self._create_radialGrid(verbose) - - elif self.ptype == 'P': - raise SyntaxError('Frequncy solver not supported for unstructured points') - #self.G = create_distArray() - - # ---------------------------------------------------------------------- - # for spatial solution, - # use the pairwise distance matrix, N x N - # where N is the number of points - elif self.stype == 'S': - - if self.ptype == 'G': - raise SyntaxError('Spatial solver not suggested for structured grids') - #self.G = self.create_points(verbose) - - elif self.ptype == 'P': - # for 1D points, add an extra axis - if (self.dim == 1) and (len(self.points.shape) == 1): - self.points = self.points[:,numpy.newaxis] - - # create pairwise distance matrix now - ltimer = Timer() - self.G = spatial.distance.cdist(self.points, self.points) - - # TODO: use a good approximation here - zval = numpy.power(10.0, -10.0) - self.G += zval * numpy.identity(self.G.shape[0]) - - if verbose: - print(' - created pairwise distance matrix:', self.G.shape, self.G.min(), self.G.max(), end = "") - ltimer.end() - - # ---------------------------------------------------------------------- - # compute the Green's function - ltimer = Timer() - if verbose: - print( ' - computing the Green\'s function:') - - if self.dim == 1: - numpy.multiply(self.G, 0.5, self.G) - - elif self.dim == 2: - numpy.log(self.G, self.G) - numpy.multiply(self.G, (0.5 / numpy.pi), self.G) - - elif self.dim == 3: - numpy.reciprocal(self.G, self.G) - numpy.multiply(self.G, (-0.25 / numpy.pi), self.G) - - if verbose: - print( self.G.shape, self.G.min(), self.G.max(), end = "") - ltimer.end() - - # ---------------------------------------------------------------------- - self.ready = True - - if verbose: - print( 'Poisson solver initialized', end ='') - gtimer.end() - - # -------------------------------------------------------------------------- - # compute the integral solution - # -------------------------------------------------------------------------- - def solve(self, f, verbose=False): - - if not self.ready: - init(verbose) - - fshape = f.shape - - gtimer = Timer() - if verbose: - print( '\nSolving Poisson Eq.') - - # ---------------------------------------------------------------------- - # regular grid - if self.ptype == 'G': - - #if (self.gdims - fshape).any(): - if self.gdims != fshape: - print( self.gdims, fshape) - raise ValueError("Shape of function should match shape of grid") - - # convolution in frequency domain - if self.stype == 'F': - p = signal.fftconvolve(f, self.G, mode='same') - numpy.multiply(p, numpy.prod(self.gdx), p) - - # convolution in spatial domain - elif self.stype == 'S': - raise SyntaxError('Spatial solver not suggested for structured grids') - #p = (self.G).dot(f.reshape(-1)).reshape(fshape) - - # ---------------------------------------------------------------------- - # unstructured points - elif self.ptype == 'P': - - if (self.points.shape[0] != fshape[0]) or (1 != len(fshape)): - raise ValueError("Shape of function should be (N,1) with N as the number of points") - - # convolution in spatial domain - # Dec 13, 2017 -- pre scaling works for unstructured grids! - if self.stype == 'S': - p = numpy.multiply(f, self.pvolumes) - numpy.dot(self.G, p, p) - #numpy.multiply(p, self.pvolumes, p) - - elif self.stype == 'F': - raise SyntaxError('Frequncy solver not supported for unstructured points') - - # ---------------------------------------------------------------------- - if verbose: - print( 'Poisson solver finished') - gtimer.end() - - return p - - # -------------------------------------------------------------------------- - # GPU-accelerated solve using CuPy - # -------------------------------------------------------------------------- - def solve_gpu(self, f, verbose=False): - """ - GPU-accelerated Poisson solver using CuPy. - Falls back to CPU if CuPy is not available. - - Args: - f: Input field (numpy array) - verbose: Print timing information - - Returns: - p: Solution (numpy array on CPU) - """ - if not CUPY_AVAILABLE: - if verbose: - print("CuPy not available, falling back to CPU solver") - return self.solve(f, verbose) - - if not self.ready: - self.prepare(verbose) - - fshape = f.shape - - gtimer = Timer() - if verbose: - print('\nSolving Poisson Eq. on GPU') - - # ---------------------------------------------------------------------- - # regular grid - if self.ptype == 'G': - if self.gdims != fshape: - print(self.gdims, fshape) - raise ValueError("Shape of function should match shape of grid") - - # convolution in frequency domain using GPU - if self.stype == 'F': - # Initialize GPU cache if needed - if self._gpu_cache is None: - if verbose: - print(' - Transferring Green\'s function to GPU...') - self._gpu_cache = { - 'G_gpu': cp.asarray(self.G, dtype=cp.float32), - 'gdx_prod': numpy.prod(self.gdx) - } - - # Transfer input to GPU - f_gpu = cp.asarray(f, dtype=cp.float32) - - # Perform FFT convolution on GPU - # CuPy's fftconvolve is GPU-accelerated - p_gpu = cp_signal.fftconvolve(f_gpu, self._gpu_cache['G_gpu'], mode='same') - - # Scale by grid spacing (on GPU) - p_gpu *= self._gpu_cache['gdx_prod'] - - # Transfer result back to CPU - p = cp.asnumpy(p_gpu) - - elif self.stype == 'S': - raise SyntaxError('Spatial solver not suggested for structured grids') - - # ---------------------------------------------------------------------- - # unstructured points - not optimized for GPU - elif self.ptype == 'P': - if verbose: - print(' - GPU solver not optimized for unstructured points, using CPU') - return self.solve(f, verbose) - - # ---------------------------------------------------------------------- - if verbose: - print('Poisson solver (GPU) finished') - gtimer.end() - - return p - - def solve_gpu_batch(self, f_list, verbose=False): - """ - Batch GPU-accelerated Poisson solver for multiple inputs. - Solves multiple Poisson equations in a single GPU call. - - Args: - f_list: List of input fields (numpy arrays) - verbose: Print timing information - - Returns: - List of solutions (numpy arrays on CPU) - """ - if not CUPY_AVAILABLE: - if verbose: - print("CuPy not available, falling back to CPU solver") - return [self.solve(f, verbose) for f in f_list] - - if not self.ready: - self.prepare(verbose) - - gtimer = Timer() - if verbose: - print(f'\nSolving {len(f_list)} Poisson Eq. on GPU (batched)') - - # Verify all inputs have same shape - fshape = f_list[0].shape - for f in f_list: - if f.shape != fshape: - raise ValueError("All input fields must have the same shape for batch processing") - - # ---------------------------------------------------------------------- - # regular grid - if self.ptype == 'G': - if self.gdims != fshape: - raise ValueError("Shape of function should match shape of grid") - - # convolution in frequency domain using GPU - if self.stype == 'F': - # Initialize GPU cache if needed - if self._gpu_cache is None: - if verbose: - print(' - Transferring Green\'s function to GPU...') - self._gpu_cache = { - 'G_gpu': cp.asarray(self.G, dtype=cp.float32), - 'gdx_prod': numpy.prod(self.gdx) - } - - results = [] - - # Process all inputs in batch - for i, f in enumerate(f_list): - # Transfer input to GPU - f_gpu = cp.asarray(f, dtype=cp.float32) - - # Perform FFT convolution on GPU - p_gpu = cp_signal.fftconvolve(f_gpu, self._gpu_cache['G_gpu'], mode='same') - - # Scale by grid spacing (on GPU) - p_gpu *= self._gpu_cache['gdx_prod'] - - # Transfer result back to CPU - results.append(cp.asnumpy(p_gpu)) - - if verbose: - print(f' - Solved {i+1}/{len(f_list)}') - - if verbose: - print('Poisson solver (GPU batch) finished') - gtimer.end() - - return results - - elif self.stype == 'S': - raise SyntaxError('Spatial solver not suggested for structured grids') - - # ---------------------------------------------------------------------- - # unstructured points - not optimized for GPU - elif self.ptype == 'P': - if verbose: - print(' - GPU solver not optimized for unstructured points, using CPU') - return [self.solve(f, verbose) for f in f_list] - - return results - - def clear_gpu_cache(self): - """Clear GPU cache to free GPU memory""" - if self._gpu_cache is not None: - self._gpu_cache = None - if CUPY_AVAILABLE: - cp.get_default_memory_pool().free_all_blocks() - - \ No newline at end of file diff --git a/src/orisys/sensor.py b/src/orisys/sensor.py index 267b3fa..906d4b3 100644 --- a/src/orisys/sensor.py +++ b/src/orisys/sensor.py @@ -17,11 +17,10 @@ from queue import Queue, Empty from .Stitcher import Stitcher from .tracker import Tracker -from .nHHD import nHHD +from .HHD import nHHD from . import configs from .get_depth import Contact3D - class Sensor: """ Optical tactile sensor interface with deformation analysis. @@ -93,21 +92,21 @@ class Sensor: self.compute_deformation_valid = False self.compute_contact_valid = False - # Data Variables - self.__img = None - self.__img_size_raw = None - self.__img_filtered = None - self.__flow = None - self.__vnormal = None - self.__vshear = None - self.__n_boxes = None - self.__fnormal = self.__fshearx = self.__fsheary = self.__ftor = 0.0 + # Data Variables (public) + self.img = None + self.img_size_raw = None + self.img_filtered = None + self.flow = None + self.vnormal = None + self.vshear = None + self.n_boxes = None + self.fnormal = self.fshearx = self.fsheary = self.ftor = 0.0 - # Contact Data Variables - self.__centroid = None - self.__centroid_found = False - self.__contour = None - self.__depth_map = None + # Contact Data Variables (public) + self.centroid = None + self.centroid_found = False + self.contour = None + self.depth_map = None # Timing self.last_log_time = time.perf_counter() @@ -182,26 +181,26 @@ class Sensor: img = np.ascontiguousarray(self.stitcher.stitch_rt_optimized()) if log_timing: timings['stitch'] = (time.perf_counter() - t0) * 1000 - self.__img_size_raw = img.shape + self.img_size_raw = img.shape # Resize img = cv2.resize(img, tuple(self._config_resolution_stitch), interpolation=cv2.INTER_CUBIC) # Temporal filtering (noise reduction) alpha = 0.7 - if self.__img_filtered is None: - self.__img_filtered = img.astype(np.float32) + if self.img_filtered is None: + self.img_filtered = img.astype(np.float32) else: - self.__img_filtered = alpha * img + (1 - alpha) * self.__img_filtered + self.img_filtered = alpha * img + (1 - alpha) * self.img_filtered - self.__img = self.__img_filtered.astype(np.uint8) + self.img = self.img_filtered.astype(np.uint8) if log_timing: timings['total'] = (time.perf_counter() - t_start) * 1000 self._print_timings("get_img", timings) if self.mirror: - self.__img = cv2.flip(self.__img, 1) + self.img = cv2.flip(self.img, 1) - return True + return self.img.copy() def compute_deformation(self): """ @@ -210,7 +209,7 @@ class Sensor: Returns: np.ndarray: Optical flow field """ - if self.__img is None: + if self.img is None: raise Exception("Warning: Get image first") self.__compute_flow() @@ -222,18 +221,17 @@ class Sensor: return True def __compute_flow(self): - if self.__img is None: + if self.img is None: raise Exception("Warning: Get image first") - flow = self.__tracker.track(self.__img) + flow = self.__tracker.track(self.img) if hasattr(flow, 'download'): flow = flow.download() - self.__flow = flow - return True + self.flow = flow def __compute_hhd(self): - flow = self.__flow + flow = self.flow if flow is None: raise Exception("Warning: Compute flow first") @@ -246,11 +244,11 @@ class Sensor: _, forces = self.__hhd.get_force_gpu_v2(flow) else: _, forces = self.__hhd.get_force(flow) - self.__fnormal, self.__fshearx, self.__fsheary, self.__ftor = forces + self.fnormal, self.fshearx, self.fsheary, self.ftor = forces # Get components - self.__vnormal = self.__hhd.get_normal() - self.__vshear = self.__hhd.get_shear() + self.vnormal = self.__hhd.get_normal() + self.vshear = self.__hhd.get_shear() # Update state self.compute_deformation_valid = True @@ -262,7 +260,7 @@ class Sensor: Compute contact information from vnormal field. Args: - vnormal: Optional vnormal field. If None, uses self.__vnormal. + vnormal: Optional vnormal field. If None, uses self.vnormal. Returns: dict: Contact result with centroid, contour, and depth_map @@ -270,7 +268,7 @@ class Sensor: multiple_foe = False # Use provided vnormal or fall back to internal state - vnormal_to_use = vnormal if vnormal is not None else self.__vnormal + vnormal_to_use = vnormal if vnormal is not None else self.vnormal if vnormal_to_use is None: raise Exception("Warning: vnormal is required for contact computation") @@ -283,11 +281,11 @@ class Sensor: else: result = self.__contact3d.compute_foe(vnormal_to_use) - self.__centroid = result['centroid'] - self.__centroid_found = result['centroid_found'] + self.centroid = result['centroid'] + self.centroid_found = result['centroid_found'] if 'contour' in result: - self.__contour = result['contour'] - self.__depth_map = result['divergence_map'] + self.contour = result['contour'] + self.depth_map = result['divergence_map'] self.compute_contact_valid = True return result @@ -304,7 +302,7 @@ class Sensor: Tuple: Requested data in order """ if any(arg in [self.info.IMG] for arg in args): - if self.__img is None: + if self.img is None: raise Exception("Warning: Get image first") if any(arg in [self.info.VRAW, self.info.VNORMAL, self.info.VSHEAR, self.info.FNORMAL, self.info.FSHEARX, self.info.FSHEARY] for arg in args): @@ -317,26 +315,26 @@ class Sensor: mapping = { # Return copies of arrays to prevent external modification of internal state - self.info.VRAW: lambda: self.__flow.copy() if self.__flow is not None else np.array([]), - self.info.VNORMAL: lambda: self.__vnormal.copy() if self.__vnormal is not None else np.array([]), - self.info.VSHEAR: lambda: self.__vshear.copy() if self.__vshear is not None else np.array([]), + self.info.VRAW: lambda: self.flow.copy() if self.flow is not None else np.array([]), + self.info.VNORMAL: lambda: self.vnormal.copy() if self.vnormal is not None else np.array([]), + self.info.VSHEAR: lambda: self.vshear.copy() if self.vshear is not None else np.array([]), # Scalar values are safe to return directly - self.info.FNORMAL: lambda: self.__fnormal, - self.info.FSHEARX: lambda: self.__fshearx, - self.info.FSHEARY: lambda: self.__fsheary, + self.info.FNORMAL: lambda: self.fnormal, + self.info.FSHEARX: lambda: self.fshearx, + self.info.FSHEARY: lambda: self.fsheary, self.info.FPS: lambda: self._fps, self.info.TIMESTAMP: lambda: time.strftime('%H:%M:%S'), # Return copies of arrays to prevent external modification - self.info.DEPTH: lambda: self.__depth_map.copy() if self.__depth_map is not None else np.array([]), + self.info.DEPTH: lambda: self.depth_map.copy() if self.depth_map is not None else np.array([]), # Centroid: can be tuple, list, or array - return a copy to prevent modification # Return None if centroid is not found or empty self.info.CENTROID: lambda: ( - np.array(self.__centroid).copy() - if self.__centroid is not None and len(self.__centroid) > 0 + np.array(self.centroid).copy() + if self.centroid is not None and len(self.centroid) > 0 else None ), - self.info.CONTOUR: lambda: self.__contour.copy() if self.__contour is not None else np.array([]), - self.info.IMG: lambda: self.__img.copy() if self.__img is not None else np.array([]), + self.info.CONTOUR: lambda: self.contour.copy() if self.contour is not None else np.array([]), + self.info.IMG: lambda: self.img.copy() if self.img is not None else np.array([]), } return tuple(mapping.get(arg, lambda: None)() for arg in args if isinstance(arg, self.info)) @@ -359,12 +357,14 @@ class Sensor: # Public API - Configuration & Calibration # ========================================================================= - def export_config(self, path="./config/exported_config.json"): - """Return current configuration dict.""" - configs.export_config(self._config, path) - if self.verbose: - print(f"Config exported to {path}") - return True + # def export_config(self, path="./config/exported_config.json"): + # """Return current configuration dict.""" + # configs.export_config(self._config, path) + # if self.verbose: + # print(f"Config exported to {path}") + # return True + + def stitch_cal(self): """ @@ -476,7 +476,7 @@ class Sensor: break v_mean0 = v_mean - self.__n_boxes = n_boxes + self.n_boxes = n_boxes self.boxes = boxes self.stitch_valid = self.isstitch and n_boxes == self._config_units @@ -609,9 +609,9 @@ class Sensor: # continue # # Put the processed image in the queue - # # get_img() updates self.__img, so we use that - # if self.__img is not None: - # self._queue_put(self._capture_queue, self.__img.copy()) + # # get_img() updates self.img, so we use that + # if self.img is not None: + # self._queue_put(self._capture_queue, self.img.copy()) # except Exception as e: # if self.verbose: print(f"Capture error: {e}") @@ -621,9 +621,9 @@ class Sensor: # try: # img = self._capture_queue.get(timeout=0.1) # # Set __img to reuse __compute_flow() - # self.__img = img + # self.img = img # self.__compute_flow() - # flow = self.__flow.copy() if self.__flow is not None else None + # flow = self.flow.copy() if self.flow is not None else None # self._queue_put(self._flow_queue, (img, flow)) # except Empty: # continue @@ -637,16 +637,16 @@ class Sensor: # img, flow = self._flow_queue.get(timeout=0.1) # # Set __flow to reuse __compute_hhd() - # self.__flow = flow + # self.flow = flow # self.__compute_hhd() # # Get results from internal variables set by __compute_hhd() # result = { # 'img': img, 'flow': flow, - # 'normal': self.__vnormal.copy() if self.__vnormal is not None else None, - # 'shear': self.__vshear.copy() if self.__vshear is not None else None, + # 'normal': self.vnormal.copy() if self.vnormal is not None else None, + # 'shear': self.vshear.copy() if self.vshear is not None else None, # # 'depth': getFilteredDensity(flow, use_cuda=self.cuda), - # 'fnormal': self.__fnormal, 'fshearx': self.__fshearx, - # 'fsheary': self.__fsheary, 'ftor': self.__ftor, + # 'fnormal': self.fnormal, 'fshearx': self.fshearx, + # 'fsheary': self.fsheary, 'ftor': self.ftor, # } # # Always include contact fields to ensure consistent result structure diff --git a/src/orisys/timer.py b/src/orisys/timer.py deleted file mode 100644 index b24f2df..0000000 --- a/src/orisys/timer.py +++ /dev/null @@ -1,59 +0,0 @@ -''' -Copyright (c) 2015, Harsh Bhatia (bhatia4@llnl.gov) -All rights reserved. - -Redistribution and use in source and binary forms, with or without -modification, are permitted provided that the following conditions are met: -1. Redistributions of source code must retain the above copyright notice, this -list of conditions and the following disclaimer. -2. Redistributions in binary form must reproduce the above copyright notice, -this list of conditions and the following disclaimer in the documentation -and/or other materials provided with the distribution. -THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND -ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED -WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE -DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE -FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL -DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER -CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, -OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -''' - -import numpy -import timeit - -class Timer(object): - - def display(self): - - tseconds = self.etime - self.stime - - print(' [[ elapsed time:', end = "") - if tseconds < numpy.power(10.,-6.): - print("%.3f micro-sec." % (tseconds*numpy.power(10,6)),end = "") - - elif tseconds < numpy.power(10.,-3.): - print("%.3f milli-sec." % (tseconds*numpy.power(10.,3.)),end = "") - - elif tseconds < 60.0: - print("%.3f sec." % (tseconds),end = "") - - else: - m = int(tseconds/ 60.0) - print("%d min. %.3f sec." % (m, (tseconds - 60*m)),end = "") - - print(']]') - - def start(self): - self.stime = timeit.default_timer() - - def end(self, print_time=True): - self.etime = timeit.default_timer() - - if print_time: - self.display() - - def __init__(self): - self.start() diff --git a/src/orisys/tracker.py b/src/orisys/tracker.py index 88dcf27..8458107 100644 --- a/src/orisys/tracker.py +++ b/src/orisys/tracker.py @@ -1,6 +1,6 @@ import cv2 import numpy as np -from .nHHD import * +# from .HHD import * # import src.utils as utils @@ -172,10 +172,10 @@ class Tracker: self.__prev_fine_flow) return self.__prev_fine_flow - def track_2dto3d(self, flow): - decompose_obj = nHHD(grid=(flow.shape[0], flow.shape[1]), spacings=(18, 18)) - flag,output = decompose_obj.get_force(flow, verbose=0) - return output + # def track_2dto3d(self, flow): + # decompose_obj = nHHD(grid=(flow.shape[0], flow.shape[1]), spacings=(18, 18)) + # flag,output = decompose_obj.get_force(flow, verbose=0) + # return output def __getForwardBackwardError(self, I0, I1, flow): if self.cuda: diff --git a/src/orisys/util.py b/src/orisys/util.py index 9df36d4..cbeb8cd 100644 --- a/src/orisys/util.py +++ b/src/orisys/util.py @@ -103,3 +103,4 @@ def draw_contact(image, contour, centroid, point_color=(0, 255, 0), region_color cv2.circle(image, (cx, cy), 10, point_color, 2) return image +