add hhd get_force gpu

example1_basic.py

@@ -13,7 +13,7 @@ def main():
 
     while True:
         sensor1.get_img() # 提取并拼接复眼图像
-        sensor1.compute_deformation() # 计算形变和分解
+        sensor1.compute_deformation(log_timing=True) # 计算形变和分解
         # 提取信息：
         # 向量场: VRAW, VNORMAL, VSHEAR: 原始，法向，切向 
         # 帧率，时间戳: FPS, TIMESTAMP

src/orisys/Stitcher.py

@@ -23,10 +23,6 @@ class Stitcher:
         self.bound=self.nominal_ov_pos - self.margin
         # for sensor 1
         self.pos_default = [(260, 56), (577, 62), (897, 69), (1214, 75), (1530, 81), (256, 391), (570, 397), (889, 403), (1208, 409), (1523, 416), (248, 724), (563, 730), (882, 738), (1201, 743), (1516, 751)]
-        # for sensor 2
-        # self.pos_default = [(195, 61), (511, 61), (831, 61), (1150, 60), (1464, 58), (195, 396), (512, 395), (830, 394), (1149, 394), (1465, 393), (197, 731), (512, 728), (831, 730), (1151, 728), (1466, 728)]
-        # for sensor 3
-        # self.pos_default = [(207, 39), (547, 40), (887, 42), (1227, 42), (1564, 44), (207, 397), (546, 399), (884, 399), (1225, 401), (1563, 403), (200, 755), (543, 757), (883, 758), (1223, 758), (1562, 760)]
  
     def stitch_rt_optimized(self, change_sv=False, s_mean=28., v_mean=180.):
         tile_list = self.tile_list

src/orisys/nHHD.py

@@ -27,6 +27,46 @@ 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):
 
@@ -148,6 +188,106 @@ class nHHD(object):
         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:
@@ -269,6 +409,374 @@ class nHHD(object):
         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

src/orisys/sensor.py

@@ -173,35 +173,87 @@ class Sensor:
         self.imgupdate = True
         return self.img
 
-    def compute_deformation(self):
+    def compute_deformation(self, log_timing=False):
         """
         Track optical flow and decompose into force field.
         
+        Args:
+            log_timing: If True, print timing information for each computation step (default: False)
+        
         Returns:
             tuple: (flow, force) where:
                 - flow: Optical flow field (2D vector field)
                 - force: Decomposed force components from nHHD， ([normal, shear_x, shear_y, torque]). Note: torque is experimental.
 
         """
+        if log_timing:
+            timings = {}
+            t_start = time.perf_counter()
+        
+        # Step 1: Optical flow tracking
+        if log_timing:
+            t0 = time.perf_counter()
         
         flow = self.tracker.track(self.img)
         
+        if log_timing:
+            timings['optical_flow'] = (time.perf_counter() - t0) * 1000
+        
+
+        
         if hasattr(flow, 'download'):
             flow = flow.download()
         
         if self.hhd is None:
             self.hhd = nHHD(grid=(flow.shape[0], flow.shape[1]), spacings=(18, 18))
 
-        self.hhd.decompose(flow)
+        if log_timing:
+            t0 = time.perf_counter()
+        # Use GPU-accelerated force computation if CUDA is enabled
+        if self.cuda and self.CUPY_AVAILABLE:
+            # Use the most optimized GPU method with custom kernel
+            flag, [self.fnormal,self.fshearx,self.fsheary,self.ftor] = self.hhd.get_force_gpu_v2(flow)
+        else:
+            flag, [self.fnormal,self.fshearx,self.fsheary,self.ftor] = self.hhd.get_force(flow)
+        # self.hhd.decompose(flow)
+        if log_timing:
+            timings['hhd_decompose'] = (time.perf_counter() - t0) * 1000
+        
         self.imgupdate = False
         self.vupdate = True
         self.flow = flow
+        
+        if log_timing:
+            t0 = time.perf_counter()
+        
         self.normal = self.hhd.get_normal()
         self.shear = self.hhd.get_shear()
-        flag, [self.fnormal,self.fshearx,self.fsheary,self.ftor] = self.hhd.get_force(flow)
+                
+        # Step 8: Compute depth/density
+        if log_timing:
+            t0 = time.perf_counter()
+        
         self.depth = getFilteredDensity(self.flow, use_cuda=self.cuda)
+        
+        if log_timing:
+            timings['get_density'] = (time.perf_counter() - t0) * 1000
+        
+        # Calculate FPS
         self.fps = 1/(time.time() - self.last_log_time)
         self.last_log_time = time.time()
+        
+        # Print timing summary
+        if log_timing:
+            timings['total'] = (time.perf_counter() - t_start) * 1000
+            print("\n[compute_deformation timing]")
+            print(f"  optical_flow:   {timings['optical_flow']:7.2f} ms")
+            print(f"  hhd_decompose:  {timings['hhd_decompose']:7.2f} ms")
+            print(f"  get_density:    {timings['get_density']:7.2f} ms")
+            print(f"  ─────────────────────────")
+            print(f"  TOTAL:          {timings['total']:7.2f} ms ({1000/timings['total']:.1f} FPS)")
+            
+            # Store timings for external access
+            self._last_timings = timings
 
         return self.flow