Merge branch 'release' into remote

.gitignore

@@ -5,3 +5,8 @@
 __pycache__/
 orisys.egg-info/
 /.cursor/
+/dist_pyd/
+/dist_single_pyd/
+/old_orisys/
+/dist_nuitka/
+/build_nuitka/

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` 方法
 
 ### 描述
 
-计算传感器变形
+按指定的信息类型（`sensor.info`枚举值）批量读取传感器状态数据，并按传入顺序返回对应值的元组。
-
+该方法通常用于在完成形变计算后，读取光流场、分量场、力、帧率等信息。
----
-## 4. `read_info` 方法
-
-### 描述
-
-获取指定类型的传感器数据，返回数量和顺序与输入参数一致。
 
 ### 输入参数
 
 * **args**: 任意数量的 `Sensor.info` 枚举，用于指定需要获取的数据类型。支持的枚举值及对应数据如下：
 
-    * VRAW: Optional[np.ndarray]          # 原始光流场, shape=(400, 400, 2)
+### 支持的枚 info keys：
-    * VNORMAL: Optional[np.ndarray]       # 法向光流场, shape=(400, 400, 2)
+向量/图像类
-    * VSHEAR: Optional[np.ndarray]            # 切向光流场, shape=(400, 400,2)
+* VRAW: 原始光流场, `np.ndarray`，shape=(H, W, 2)
-    * DEPTH: Optional[np.ndarray]        # 深度场, shape=(400, 400,1)
+* 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]            # 横向切向力分量
+* FNORMAL: 法向力，`float`。
-    * FSHEARY: Optional[np.float64]        # 纵向切向力分量
+* 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
+### 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：若未指定或未找到标定文件，系统会在需要时创建新的拼接标定数据并保存到默认路径。
+拼接结果异常通常与标定数据、拼接配置或输入图像质量有关，建议重新生成拼接标定并检查配置参数。

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,22 +158,23 @@ 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}")
+            configs_dst.mkdir(parents=True, exist_ok=True)
-            shutil.copytree(configs_src, configs_dst, dirs_exist_ok=True)
+        
-        else:
+        # Only copy JSON files, not any .py files
-            # 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)")
+            print(f"Copying JSON files from configs directory ({len(json_files_src)} files)")
-                shutil.copytree(configs_src, configs_dst, dirs_exist_ok=True)
+            for json_file in json_files_src:
+                shutil.copy(json_file, configs_dst / json_file.name)
         else:
             print(f"Configs directory already present with {len(json_files_dst)} JSON files")
     
@@ -157,12 +184,46 @@ 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):
     """Create a modified setup.py for the obfuscated package."""
@@ -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__':

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()
+

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()
+

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' 重置光流追踪器")
 

src/orisys/MultiSensor.py

@@ -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()
-

src/orisys/StructuredGrid.py

@@ -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()
-
-    # --------------------------------------------------------------------------
-    # --------------------------------------------------------------------------

src/orisys/UnstructuredGrid.py

@@ -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
-
-    # --------------------------------------------------------------------------
-    # --------------------------------------------------------------------------

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'
+]

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.
@@ -97,25 +170,3 @@ def _format_json_value(value, indent_level=0, indent_str="    "):
     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
-    

src/orisys/getDensity.py

@@ -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()

src/orisys/nHHD.py

@@ -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)
-# ---------------------------------------------------------------------------

src/orisys/nPoisson.py

@@ -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()
-    
-    

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
+        # Data Variables (public)
-        self.__img = None
+        self.img = None
-        self.__img_size_raw = None
+        self.img_size_raw = None
-        self.__img_filtered = None
+        self.img_filtered = None
-        self.__flow = None
+        self.flow = None
-        self.__vnormal = None
+        self.vnormal = None
-        self.__vshear = None
+        self.vshear = None
-        self.__n_boxes = None
+        self.n_boxes = None
-        self.__fnormal = self.__fshearx = self.__fsheary = self.__ftor = 0.0
+        self.fnormal = self.fshearx = self.fsheary = self.ftor = 0.0
 
-        # Contact Data Variables
+        # Contact Data Variables (public)
-        self.__centroid = None
+        self.centroid = None
-        self.__centroid_found = False
+        self.centroid_found = False
-        self.__contour = None
+        self.contour = None
-        self.__depth_map = 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:
+        if self.img_filtered is None:
-            self.__img_filtered = img.astype(np.float32)
+            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
+        self.flow = flow
-        return True
     
     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.vnormal = self.__hhd.get_normal()
-        self.__vshear = self.__hhd.get_shear()
+        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 = result['centroid']
-        self.__centroid_found = result['centroid_found']
+        self.centroid_found = result['centroid_found']
         if 'contour' in result:
-            self.__contour = result['contour']
+            self.contour = result['contour']
-        self.__depth_map = result['divergence_map']
+        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.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.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.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.FNORMAL: lambda: self.fnormal,
-            self.info.FSHEARX: lambda: self.__fshearx,
+            self.info.FSHEARX: lambda: self.fshearx,
-            self.info.FSHEARY: lambda: self.__fsheary,
+            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() 
+                np.array(self.centroid).copy() 
-                if self.__centroid is not None and len(self.__centroid) > 0 
+                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.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.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"):
+    # def export_config(self, path="./config/exported_config.json"):
-        """Return current configuration dict."""
+    #     """Return current configuration dict."""
-        configs.export_config(self._config, path)
+    #     configs.export_config(self._config, path)
-        if self.verbose:
+    #     if self.verbose:
-            print(f"Config exported to {path}")
+    #         print(f"Config exported to {path}")
-        return True
+    #     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
+    #             # get_img() updates self.img, so we use that
-    #             if self.__img is not None:
+    #             if self.img is not None:
-    #                 self._queue_put(self._capture_queue, self.__img.copy())
+    #                 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,
+    #                 'normal': self.vnormal.copy() if self.vnormal is not None else None,
-    #                 'shear': self.__vshear.copy() if self.__vshear 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,
+    #                 'fnormal': self.fnormal, 'fshearx': self.fshearx,
-    #                 'fsheary': self.__fsheary, 'ftor': self.__ftor,
+    #                 'fsheary': self.fsheary, 'ftor': self.ftor,
     #             }
                 
     #             # Always include contact fields to ensure consistent result structure

src/orisys/timer.py

@@ -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()

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):
+    # def track_2dto3d(self, flow):
-        decompose_obj = nHHD(grid=(flow.shape[0], flow.shape[1]), spacings=(18, 18))
+    #     decompose_obj = nHHD(grid=(flow.shape[0], flow.shape[1]), spacings=(18, 18))
-        flag,output = decompose_obj.get_force(flow, verbose=0)
+    #     flag,output = decompose_obj.get_force(flow, verbose=0)
-        return output
+    #     return output
         
     def __getForwardBackwardError(self, I0, I1, flow):
         if self.cuda:

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
+