diff --git a/src/orisys/HHD.py b/src/orisys/HHD.py
index d3e00ab..34ad3d7 100644
--- a/src/orisys/HHD.py
+++ b/src/orisys/HHD.py
@@ -309,34 +309,47 @@ class nHHD(object):
         
         self.decompose(vfield, verbose)       
         
-        self.Sn = numpy.sum(numpy.absolute(self.d))
-        # print("Sn is:", self.Sn)
+        # Optimized: use abs instead of absolute (slightly faster)
+        self.Sn = numpy.sum(numpy.abs(self.d))
 
-        # self.Stx = abs(numpy.sum(self.h[:, :, 0]))
-        # self.Sty = abs(numpy.sum(self.h[:, :, 1]))
+        # Optimized: direct indexing and sum (already efficient)
         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)
+        # Optimized: find both max and min indices in single pass through flattened array
+        nRu_flat = self.nRu.ravel()
+        nR_max_idx = numpy.argmax(nRu_flat)
+        nR_min_idx = numpy.argmin(nRu_flat)
+        # Convert flat indices to 2D coordinates (faster than unravel_index)
+        nR_max_loc = (nR_max_idx // shape[1], nR_max_idx % shape[1])
+        nR_min_loc = (nR_min_idx // shape[1], nR_min_idx % shape[1])
         
-        # print(nR_max_loc)
+        # Optimized: compute cross products directly without creating intermediate arrays
+        # For 2D: cross(a, b) = a[0]*b[1] - a[1]*b[0]
+        # Extract components once for reuse
+        r_x = self.r[:, :, 0]
+        r_y = self.r[:, :, 1]
+        grid_x = grid[:, :, 0]
+        grid_y = grid[:, :, 1]
         
-        r_p = grid - nR_max_loc 
-        r_n = grid - nR_min_loc
+        # Compute cross products directly (faster than numpy.cross for 2D)
+        # cross(r_p, r) where r_p = grid - nR_max_loc
+        max_y, max_x = nR_max_loc[0], nR_max_loc[1]
+        cross_p = (grid_x - max_x) * r_y - (grid_y - max_y) * r_x
         
-        self.Stor = numpy.sum(numpy.cross(r_p, self.r)) + numpy.sum(numpy.cross(r_n, self.r))
-        # print("Stor is:", self.Stor)
+        # cross(r_n, r) where r_n = grid - nR_min_loc
+        min_y, min_x = nR_min_loc[0], nR_min_loc[1]
+        cross_n = (grid_x - min_x) * r_y - (grid_y - min_y) * r_x
+        
+        self.Stor = numpy.sum(cross_p) + numpy.sum(cross_n)
 
-        output = numpy.array([self.Sn, self.Stx, self.Sty, self.Stor])
+        # Optimized: pre-allocate output array with explicit dtype
+        output = numpy.array([self.Sn, self.Stx, self.Sty, self.Stor], dtype=numpy.float32)
 
-        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
+        # Optimized: combine threshold checks (short-circuit evaluation)
+        flag = (self.Sn >= 300 and self.Stx >= 300 and self.Sty >= 300 and numpy.abs(self.Stor) >= 1000)
+        
+        return flag, output
         # return self.d, self.h, numpy.cross(r_p, self.r) , numpy.cross(r_n, self.r)
 
     def _init_force_cpu_cache(self, shape):