atrstm

commit 477fb701f096b706c222a18958146814be6d5a8b
parent 84945198e83707fb23388d4858849a97f4346fdf
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Wed, 27 Jan 2021 14:55:59 +0100

Upd the auto-computation of the "grid definition"

Diffstat:
M
src/atrstm.c
 | 
34
+++++++++++++++++-----------------

1 file changed, 17 insertions(+), 17 deletions(-)
diff --git a/src/atrstm.c b/src/atrstm.c
@@ -64,6 +64,17 @@ check_args(const struct atrstm_args* args)
   return 1;
 }
 
+static FINLINE unsigned
+round_pow2(const unsigned val)
+{
+  const unsigned next_pow2 = (unsigned)round_up_pow2(val);
+  if(next_pow2 - val <= next_pow2/4) {
+    return next_pow2;
+  } else {
+    return next_pow2/2;
+  }
+}
+
 static void
 compute_default_grid_definition
   (const struct atrstm* atrstm,
@@ -77,7 +88,6 @@ compute_default_grid_definition
   double vxsz;
   int iaxis_max;
   int iaxis_remain[2];
-  size_t next_pow2[2];
   int i;
   ASSERT(atrstm && def);
 
@@ -98,26 +108,16 @@ compute_default_grid_definition
 
   /* Fix the definition to along the maximum axis and compute the voxel size
    * along this axis */
-  def[iaxis_max] = (unsigned)round_up_pow2(definition_hint);
+  def[iaxis_max] = round_pow2(definition_hint);
   vxsz = sz[iaxis_max] / def[iaxis_max];
 
-  /* Compute the definition along the 2 remaining axes. First, copute these
-   * definitions by assuming that the voxels ar cubics [...] */
+  /* Compute the definition along the 2 remaining axes. First, compute them by
+   * assuming that the voxels are cubics and then round these definitions to
+   * their nearest power of two. */
   def[iaxis_remain[0]] = (unsigned)lround(sz[iaxis_remain[0]]/vxsz);
   def[iaxis_remain[1]] = (unsigned)lround(sz[iaxis_remain[1]]/vxsz);
-
-  /* [...] Next compute the power of two directly greatest that the computed
-   * definition [...] */
-  next_pow2[0] = round_up_pow2(def[iaxis_remain[0]]);
-  next_pow2[1] = round_up_pow2(def[iaxis_remain[1]]);
-
-  /* [...] Finally round the definition to its nearest power of two. */
-  def[iaxis_remain[0]] =
-    next_pow2[0] - def[iaxis_remain[0]] < next_pow2[0]/4
-    ? (unsigned)next_pow2[0] : (unsigned)next_pow2[0]/2;
-  def[iaxis_remain[1]] =
-    next_pow2[1] - def[iaxis_remain[1]] < next_pow2[1]/4
-    ? (unsigned)next_pow2[1] : (unsigned)next_pow2[1]/2;
+  def[iaxis_remain[0]] = round_pow2(def[iaxis_remain[0]]);
+  def[iaxis_remain[1]] = round_pow2(def[iaxis_remain[1]]);
 }
 
 static void