diff --git a/include/openmc/lattice.h b/include/openmc/lattice.h
index ba7cbdf5e..9c8b8af68 100644
--- a/include/openmc/lattice.h
+++ b/include/openmc/lattice.h
@@ -104,7 +104,7 @@ public:
   //! \brief Find the lattice tile indices for a given point.
   //! \param r A 3D Cartesian coordinate.
   //! \return An array containing the indices of a lattice tile.
-  virtual std::array<int, 3> get_indices(Position r) const = 0;
+  virtual std::array<int, 3> get_indices(Position r, Direction u) const = 0;
 
   //! \brief Get coordinates local to a lattice tile.
   //! \param r A 3D Cartesian coordinate.
@@ -212,7 +212,7 @@ public:
   std::pair<double, std::array<int, 3>>
   distance(Position r, Direction u, const std::array<int, 3>& i_xyz) const;
 
-  std::array<int, 3> get_indices(Position r) const;
+  std::array<int, 3> get_indices(Position r, Direction u) const;
 
   Position
   get_local_position(Position r, const std::array<int, 3> i_xyz) const;
@@ -252,7 +252,7 @@ public:
   std::pair<double, std::array<int, 3>>
   distance(Position r, Direction u, const std::array<int, 3>& i_xyz) const;
 
-  std::array<int, 3> get_indices(Position r) const;
+  std::array<int, 3> get_indices(Position r, Direction u) const;
 
   Position
   get_local_position(Position r, const std::array<int, 3> i_xyz) const;
diff --git a/src/geometry.cpp b/src/geometry.cpp
index feaf6c00f..3c7799c27 100644
--- a/src/geometry.cpp
+++ b/src/geometry.cpp
@@ -215,8 +215,7 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
       // Determine lattice indices.
       Position r {p->coord_[p->n_coord_-1].xyz};
       Direction u {p->coord_[p->n_coord_-1].uvw};
-      r += TINY_BIT * u;
-      auto i_xyz = lat.get_indices(r);
+      auto i_xyz = lat.get_indices(r, u);
 
       // Store lower level coordinates.
       r = lat.get_local_position(p->coord_[p->n_coord_-1].xyz, i_xyz);
diff --git a/src/lattice.cpp b/src/lattice.cpp
index 780300b28..6e0f2331f 100644
--- a/src/lattice.cpp
+++ b/src/lattice.cpp
@@ -285,15 +285,38 @@ const
 //==============================================================================
 
 std::array<int, 3>
-RectLattice::get_indices(Position r) const
+RectLattice::get_indices(Position r, Direction u) const
 {
-  int ix {static_cast<int>(std::ceil((r.x - lower_left_.x) / pitch_.x))-1};
-  int iy {static_cast<int>(std::ceil((r.y - lower_left_.y) / pitch_.y))-1};
-  int iz;
-  if (is_3d_) {
-    iz = static_cast<int>(std::ceil((r.z - lower_left_.z) / pitch_.z))-1;
+  // Determine x index, accounting for coincidence
+  double ix_ {(r.x - lower_left_.x) / pitch_.x};
+  long ix_close {std::lround(ix_)};
+  int ix;
+  if (std::abs(ix_ - ix_close) < FP_COINCIDENT) {
+    ix = (u.x > 0) ? ix_close : ix_close - 1;
   } else {
-    iz = 0;
+    ix = std::floor(ix_);
+  }
+
+  // Determine y index, accounting for coincidence
+  double iy_ {(r.y - lower_left_.y) / pitch_.y};
+  long iy_close {std::lround(iy_)};
+  int iy;
+  if (std::abs(iy_ - iy_close) < FP_COINCIDENT) {
+    iy = (u.y > 0) ? iy_close : iy_close - 1;
+  } else {
+    iy = std::floor(iy_);
+  }
+
+  // Determine z index, accounting for coincidence
+  int iz = 0;
+  if (is_3d_) {
+    double iz_ {(r.z - lower_left_.z) / pitch_.z};
+    long iz_close {std::lround(iz_)};
+    if (std::abs(iz_ - iz_close) < FP_COINCIDENT) {
+      iz = (u.z > 0) ? iz_close : iz_close - 1;
+    } else {
+      iz = std::floor(iz_);
+    }
   }
   return {ix, iy, iz};
 }
@@ -685,8 +708,13 @@ const
 //==============================================================================
 
 std::array<int, 3>
-HexLattice::get_indices(Position r) const
+HexLattice::get_indices(Position r, Direction u) const
 {
+  // The implementation for HexLattice currently doesn't use direction
+  // information. As a result, we move the position slightly forward to
+  // determine what lattice index the particle is most likely to be in.
+  r += TINY_BIT * u;
+
   // Offset the xyz by the lattice center.
   Position r_o {r.x - center_.x, r.y - center_.y, r.z};
   if (is_3d_) {r_o.z -= center_.z;}
@@ -694,7 +722,7 @@ HexLattice::get_indices(Position r) const
   // Index the z direction.
   std::array<int, 3> out;
   if (is_3d_) {
-    out[2] = static_cast<int>(std::ceil(r_o.z / pitch_[1] + 0.5 * n_axial_))-1;
+    out[2] = std::floor(r_o.z / pitch_[1] + 0.5 * n_axial_);
   } else {
     out[2] = 0;
   }
@@ -702,9 +730,8 @@ HexLattice::get_indices(Position r) const
   // Convert coordinates into skewed bases.  The (x, alpha) basis is used to
   // find the index of the global coordinates to within 4 cells.
   double alpha = r_o.y - r_o.x / std::sqrt(3.0);
-  out[0] = static_cast<int>(std::floor(r_o.x
-                                       / (0.5*std::sqrt(3.0) * pitch_[0])));
-  out[1] = static_cast<int>(std::floor(alpha / pitch_[0]));
+  out[0] = std::floor(r_o.x / (0.5*std::sqrt(3.0) * pitch_[0]));
+  out[1] = std::floor(alpha / pitch_[0]);
 
   // Add offset to indices (the center cell is (i_x, i_alpha) = (0, 0) but
   // the array is offset so that the indices never go below 0).
@@ -737,12 +764,12 @@ HexLattice::get_indices(Position r) const
   // Select the minimum squared distance which corresponds to the cell the
   // coordinates are in.
   if (k_min == 2) {
-    out[0] += 1;
+    ++out[0];
   } else if (k_min == 3) {
-    out[1] += 1;
+    ++out[1];
   } else if (k_min == 4) {
-    out[0] += 1;
-    out[1] += 1;
+    ++out[0];
+    ++out[1];
   }
 
   return out;