# Intersection of region with triangular lattice in 2D

Is there a function that finds all points that lie in specified region and also are part of triangular lattice defined by {Sqrt[3] i, 3 (2 j + Mod[i, 2])} for integer i, j?

reg = Region@
Polygon@{{0, 0}, {4 Sqrt[3], 0}, {8 Sqrt[3], 12}, {2 Sqrt[3],
30}, {-2 Sqrt[3], 30}, {-6 Sqrt[3], 18}};
poi = Point[{{0, 0}, {2 Sqrt[3], 0}, {4 Sqrt[3], 0}, {-Sqrt[3],
3}, {Sqrt[3], 3}, {3 Sqrt[3], 3}, {5 Sqrt[3], 3}, {-2 Sqrt[3],
6}, {0, 6}, {2 Sqrt[3], 6}, {4 Sqrt[3], 6}, {6 Sqrt[3],
6}, {-3 Sqrt[3], 9}, {-Sqrt[3], 9}, {Sqrt[3], 9}, {3 Sqrt[3],
9}, {5 Sqrt[3], 9}, {7 Sqrt[3], 9}, {-4 Sqrt[3],
12}, {-2 Sqrt[3], 12}, {0, 12}, {2 Sqrt[3], 12}, {4 Sqrt[3],
12}, {6 Sqrt[3], 12}, {8 Sqrt[3], 12}, {-5 Sqrt[3],
15}, {-3 Sqrt[3], 15}, {-Sqrt[3], 15}, {Sqrt[3], 15}, {3 Sqrt[3],
15}, {5 Sqrt[3], 15}, {7 Sqrt[3], 15}, {-6 Sqrt[3],
18}, {-4 Sqrt[3], 18}, {-2 Sqrt[3], 18}, {0, 18}, {2 Sqrt[3],
18}, {4 Sqrt[3], 18}, {6 Sqrt[3], 18}, {-5 Sqrt[3],
21}, {-3 Sqrt[3], 21}, {-Sqrt[3], 21}, {Sqrt[3], 21}, {3 Sqrt[3],
21}, {5 Sqrt[3], 21}, {-4 Sqrt[3], 24}, {-2 Sqrt[3], 24}, {0,
24}, {2 Sqrt[3], 24}, {4 Sqrt[3], 24}, {-3 Sqrt[3],
27}, {-Sqrt[3], 27}, {Sqrt[3], 27}, {3 Sqrt[3], 27}, {-2 Sqrt[3],
30}, {0, 30}, {2 Sqrt[3], 30}}];
Show[reg, Graphics[{Black, poi}]]
{Sqrt[3] i, 3 (2 j + Mod[i, 2])}


In other words you are given region reg and lattice defined by {Sqrt[3] i, 3 (2 j + Mod[i, 2])} for integer i, j. How to find all those points which I gave in poi?

I can do it by cycling over all possible i, j and using RegionMember but built in function would be better.

Update:

In fact it can be done by the following code, but it is very sloooow. Cycling over i, j is much faster.

Solve[RegionMember[reg, {Sqrt[3] i, 3 (2 j + Mod[i, 2])}], Integers];
{Sqrt[3] i, 3 (2 j + Mod[i, 2])} /. %;
Sort[Identity @@ poi] == Sort[%]

*( True )*


The Solve had probably problems with {Sqrt[3] i, 3 (2 j + Mod[i, 2])} but if I use {Sqrt[3] (2 x - y), 3 y} (which defines the same triangular lattice) the computation is much faster.
Solve[RegionMember[reg, {Sqrt[3] (2 x - y), 3 y}], Integers];