5#ifndef GEMMI_ASUMASK_HPP_
6#define GEMMI_ASUMASK_HPP_
16 static constexpr int denom = 24;
26 std::string
str()
const {
28 for (
int i = 0;
i < 3; ++
i) {
33 s +=
incl[
i] ?
"<=" :
"<";
34 impl::append_fraction(s, impl::get_op_fraction(
size[
i]));
56 fail(
"grid is not fully setup");
59 auto iceil = [](
double x) {
return int(x) + 1; };
69 fail(
"Missing space group");
71 using Point = std::array<int, 3>;
85 const std::vector<GridOp>
ops = grid.get_scaled_ops_except_id();
129 auto is_in = [&](
const Point&
p) {
139 grid.index_n_ref(
t[0],
t[1],
t[2]);
148 for (;
it != end; ++
it)
158 for (
int w = 0; w <
w_lim; ++w)
159 for (
int v = 0; v <
v_lim; ++v)
160 for (
int u = 0; u <
u_lim; ++u) {
161 size_t idx = grid.index_q(u, v, w);
162 if (grid.data[idx] == 0) {
165 Point
t = op.apply(u, v, w);
174 for (
size_t n = 0; n != grid.data.size(); ++n)
175 if (grid.data[n] == 0) {
176 auto p = grid.index_to_point(n);
177 printf(
"[debug1] %d %d %d is missing\n",
p.u,
p.v,
p.w);
181 printf(
"[debug2] checkpoints failed\n");
183 if (std::find(grid.data.begin(), grid.data.end(), 0) != grid.data.end())
202 return a.volume < b.volume ||
203 (a.volume == b.volume && a.size[0] + a.size[1] + a.size[2] <
204 b.size[0] + b.size[1] + b.size[2]);
208 for (
int i = 0;
i < 3; ++
i) {
221template<
typename T,
typename V=std::
int8_t>
struct MaskedGrid {
254template<
typename V=std::
int8_t>
259 for (
int w = 0; w < end[2]; ++w)
260 for (
int v = 0; v < end[1]; ++v)
261 for (
int u = 0; u < end[0]; ++u) {
262 size_t idx = grid.
index_q(u, v, w);
263 if (mask[idx] == 2) {
266 std::array<int, 3>
t = op.apply(u, v, w);
274 if (std::find(mask.begin(), mask.end(), 2) != mask.end())
275 fail(
"get_asu_mask(): internal error");
289inline std::pair<int, int> trim_false_values(
const std::vector<bool>& vec) {
290 const int n = (int) vec.size();
292 std::pair<int, int> span{n, n};
295 if (!vec[0] || !vec[n-1]) {
297 while (span.first != 0 && !vec[span.first-1])
303 while (span.second != n && !vec[span.second])
305 max_trim = span.second + (n - span.first);
307 for (
int start = 0; ;) {
316 while (end != n && !vec[end])
318 if (end - start > max_trim) {
319 max_trim = end - start;
332 grid.check_not_empty();
333 std::vector<bool> nonzero[3];
334 nonzero[0].resize(grid.nu,
false);
335 nonzero[1].resize(grid.nv,
false);
336 nonzero[2].resize(grid.nw,
false);
338 for (
int w = 0; w != grid.nw; ++w)
339 for (
int v = 0; v != grid.nv; ++v)
340 for (
int u = 0; u != grid.nu; ++u, ++idx) {
341 T val = grid.data[idx];
342 if (!(impl::is_nan(val) || val == 0)) {
343 nonzero[0][u] =
true;
344 nonzero[1][v] =
true;
345 nonzero[2][w] =
true;
349 for (
int i = 0;
i < 3; ++
i) {
350 auto span = impl::trim_false_values(nonzero[
i]);
351 double inv_n = 1.0 / nonzero[
i].size();
352 box.minimum.at(
i) = (span.second - 0.5) *
inv_n -
int(span.second >= span.first);
353 box.maximum.at(
i) = (span.first - 0.5) *
inv_n;
MaskedGrid< T > masked_asu(Grid< T > &grid)
AsuBrick find_asu_brick(const SpaceGroup *sg)
Box< Fractional > get_nonzero_extent(const GridBase< T > &grid)
void fail(const std::string &msg)
std::vector< V > get_asu_mask(const GridMeta &grid)
std::array< int, 3 > size
Box< Fractional > get_extent() const
Fractional get_upper_limit() const
std::array< bool, 3 > incl
static constexpr int denom
std::array< int, 3 > uvw_end(const GridMeta &meta) const
AsuBrick(int a, int b, int c)
grid coordinates (modulo size) and a pointer to value
bool operator!=(const iterator &o) const
iterator(MaskedGrid &parent_, size_t index_)
GridBase< T >::Point operator*()
bool operator==(const iterator &o) const