Skip to content
GitLab
Menu
Projects
Groups
Snippets
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
Simon Perkins
ducc
Commits
8e9b2689
Commit
8e9b2689
authored
Jan 27, 2020
by
Martin Reinecke
Browse files
datatype patrol
parent
09ce0bd3
Changes
5
Hide whitespace changes
Inline
Side-by-side
Healpix_cxx/healpix_base.cc
View file @
8e9b2689
...
...
@@ -100,9 +100,9 @@ namespace {
2: pixel center is inside the shape, but maybe not the complete pixel
3: pixel lies completely inside the shape */
template
<
typename
I
,
typename
I2
>
inline
void
check_pixel
(
int
o
,
in
t
order_
,
in
t
omax
,
in
t
zone
,
rangeset
<
I2
>
&
pixset
,
I
pix
,
vector
<
pair
<
I
,
in
t
>
>
&
stk
,
bool
inclusive
,
in
t
&
stacktop
)
template
<
typename
I
,
typename
I2
>
inline
void
check_pixel
(
size_t
o
,
size_
t
order_
,
size_
t
omax
,
size_
t
zone
,
rangeset
<
I2
>
&
pixset
,
I
pix
,
vector
<
pair
<
I
,
size_
t
>
>
&
stk
,
bool
inclusive
,
size_
t
&
stacktop
)
{
if
(
zone
==
0
)
return
;
...
...
@@ -114,7 +114,7 @@ template<typename I, typename I2> inline void check_pixel (int o, int order_,
pixset
.
append
(
pix
<<
sdist
,(
pix
+
1
)
<<
sdist
);
// output all subpixels
}
else
// (1<=zone<=2)
for
(
in
t
i
=
0
;
i
<
4
;
++
i
)
for
(
size_
t
i
=
0
;
i
<
4
;
++
i
)
stk
.
push_back
(
make_pair
(
4
*
pix
+
3
-
i
,
o
+
1
));
// add children
}
else
if
(
o
>
order_
)
// this implies that inclusive==true
...
...
@@ -145,7 +145,7 @@ template<typename I, typename I2> inline void check_pixel (int o, int order_,
if
(
order_
<
omax
)
// check sublevels
{
stacktop
=
stk
.
size
();
// remember current stack position
for
(
in
t
i
=
0
;
i
<
4
;
++
i
)
// add children in reverse order
for
(
size_
t
i
=
0
;
i
<
4
;
++
i
)
// add children in reverse order
stk
.
push_back
(
make_pair
(
4
*
pix
+
3
-
i
,
o
+
1
));
}
else
// at resolution limit
...
...
@@ -326,12 +326,12 @@ template<typename I> template<typename I2>
crpdr
[
o
]
=
(
radius
+
dr
>
pi
)
?
-
1.
:
cos
(
radius
+
dr
);
crmdr
[
o
]
=
(
radius
-
dr
<
0.
)
?
1.
:
cos
(
radius
-
dr
);
}
vector
<
pair
<
I
,
in
t
>
>
stk
;
// stack for pixel numbers and their orders
vector
<
pair
<
I
,
size_
t
>
>
stk
;
// stack for pixel numbers and their orders
stk
.
reserve
(
12
+
3
*
omax
);
// reserve maximum size to avoid reallocation
for
(
int
i
=
0
;
i
<
12
;
++
i
)
// insert the 12 base pixels in reverse order
stk
.
push_back
(
make_pair
(
I
(
11
-
i
),
0
));
in
t
stacktop
=
0
;
// a place to save a stack position
size_
t
stacktop
=
0
;
// a place to save a stack position
while
(
!
stk
.
empty
())
// as long as there are pixels on the stack
{
...
...
@@ -347,7 +347,7 @@ template<typename I> template<typename I2>
if
(
cangdist
>
crpdr
[
o
])
{
in
t
zone
=
(
cangdist
<
cosrad
)
?
1
:
((
cangdist
<=
crmdr
[
o
])
?
2
:
3
);
size_
t
zone
=
(
cangdist
<
cosrad
)
?
1
:
((
cangdist
<=
crmdr
[
o
])
?
2
:
3
);
check_pixel
(
o
,
order_
,
omax
,
zone
,
pixset
,
pix
,
stk
,
inclusive
,
stacktop
);
...
...
@@ -477,7 +477,7 @@ template<typename I> template<typename I2>
}
else
// scheme_ == NEST
{
in
t
oplus
=
0
;
size_
t
oplus
=
0
;
if
(
inclusive
)
{
MR_assert
((
I
(
1
)
<<
(
order_max
-
order_
))
>=
fact
,
...
...
@@ -486,13 +486,13 @@ template<typename I> template<typename I2>
"oversampling factor must be a power of 2"
);
oplus
=
ilog2
(
fact
);
}
in
t
omax
=
order_
+
oplus
;
// the order up to which we test
size_
t
omax
=
size_t
(
order_
)
+
oplus
;
// the order up to which we test
// TODO: ignore all disks with radius>=pi
vector
<
T_Healpix_Base
<
I
>
>
base
(
omax
+
1
);
mav
<
double
,
3
>
crlimit
({
size_t
(
omax
)
+
1
,
nv
,
3
});
for
(
in
t
o
=
0
;
o
<=
omax
;
++
o
)
// prepare data at the required orders
for
(
size_
t
o
=
0
;
o
<=
omax
;
++
o
)
// prepare data at the required orders
{
base
[
o
].
Set
(
o
,
NEST
);
double
dr
=
base
[
o
].
max_pixrad
();
// safety distance
...
...
@@ -504,18 +504,18 @@ template<typename I> template<typename I2>
}
}
vector
<
pair
<
I
,
in
t
>
>
stk
;
// stack for pixel numbers and their orders
vector
<
pair
<
I
,
size_
t
>
>
stk
;
// stack for pixel numbers and their orders
stk
.
reserve
(
12
+
3
*
omax
);
// reserve maximum size to avoid reallocation
for
(
int
i
=
0
;
i
<
12
;
++
i
)
// insert the 12 base pixels in reverse order
stk
.
push_back
(
make_pair
(
I
(
11
-
i
),
0
));
in
t
stacktop
=
0
;
// a place to save a stack position
size_
t
stacktop
=
0
;
// a place to save a stack position
while
(
!
stk
.
empty
())
// as long as there are pixels on the stack
{
// pop current pixel number and order from the stack
I
pix
=
stk
.
back
().
first
;
in
t
o
=
stk
.
back
().
second
;
size_
t
o
=
stk
.
back
().
second
;
stk
.
pop_back
();
vec3
pv
(
base
[
o
].
pix2vec
(
pix
));
...
...
@@ -550,14 +550,14 @@ template<typename I> void T_Healpix_Base<I>::query_multidisc_general
}
else
// scheme_ == NEST
{
in
t
oplus
=
inclusive
?
2
:
0
;
in
t
omax
=
min
<
in
t
>
(
order_max
,
order_
+
oplus
);
// the order up to which we test
size_
t
oplus
=
inclusive
?
2
:
0
;
size_
t
omax
=
min
<
size_
t
>
(
order_max
,
size_t
(
order_
)
+
oplus
);
// the order up to which we test
// TODO: ignore all disks with radius>=pi
vector
<
T_Healpix_Base
<
I
>
>
base
(
omax
+
1
);
mav
<
double
,
3
>
crlimit
({
size_t
(
omax
+
1
),
nv
,
3
});
for
(
in
t
o
=
0
;
o
<=
omax
;
++
o
)
// prepare data at the required orders
for
(
size_
t
o
=
0
;
o
<=
omax
;
++
o
)
// prepare data at the required orders
{
base
[
o
].
Set
(
o
,
NEST
);
double
dr
=
base
[
o
].
max_pixrad
();
// safety distance
...
...
@@ -569,12 +569,12 @@ template<typename I> void T_Healpix_Base<I>::query_multidisc_general
}
}
vector
<
pair
<
I
,
in
t
>
>
stk
;
// stack for pixel numbers and their orders
vector
<
pair
<
I
,
size_
t
>
>
stk
;
// stack for pixel numbers and their orders
stk
.
reserve
(
12
+
3
*
omax
);
// reserve maximum size to avoid reallocation
for
(
int
i
=
0
;
i
<
12
;
++
i
)
// insert the 12 base pixels in reverse order
stk
.
push_back
(
make_pair
(
I
(
11
-
i
),
0
));
in
t
stacktop
=
0
;
// a place to save a stack position
size_
t
stacktop
=
0
;
// a place to save a stack position
vector
<
size_t
>
zone
(
nv
);
vector
<
size_t
>
zstk
;
zstk
.
reserve
(
cmds
.
size
());
...
...
@@ -583,7 +583,7 @@ template<typename I> void T_Healpix_Base<I>::query_multidisc_general
{
// pop current pixel number and order from the stack
I
pix
=
stk
.
back
().
first
;
in
t
o
=
stk
.
back
().
second
;
size_
t
o
=
stk
.
back
().
second
;
stk
.
pop_back
();
vec3
pv
(
base
[
o
].
pix2vec
(
pix
));
...
...
@@ -1161,7 +1161,7 @@ template<typename I> void T_Healpix_Base<I>::neighbors (I pix,
if
((
ix
>
0
)
&&
(
ix
<
nsm1
)
&&
(
iy
>
0
)
&&
(
iy
<
nsm1
))
{
if
(
scheme_
==
RING
)
for
(
in
t
m
=
0
;
m
<
8
;
++
m
)
for
(
size_
t
m
=
0
;
m
<
8
;
++
m
)
result
[
m
]
=
xyf2ring
(
ix
+
nb_xoffset
[
m
],
iy
+
nb_yoffset
[
m
],
face_num
);
else
{
...
...
@@ -1178,7 +1178,7 @@ template<typename I> void T_Healpix_Base<I>::neighbors (I pix,
}
else
{
for
(
in
t
i
=
0
;
i
<
8
;
++
i
)
for
(
size_
t
i
=
0
;
i
<
8
;
++
i
)
{
int
x
=
ix
+
nb_xoffset
[
i
],
y
=
iy
+
nb_yoffset
[
i
];
int
nbnum
=
4
;
...
...
Healpix_cxx/healpix_tables.cc
View file @
8e9b2689
...
...
@@ -104,9 +104,9 @@ const int Healpix_Tables::nb_swaparray[][3] =
{
6
,
0
,
0
},
// NW
{
3
,
0
,
0
}
};
// N
const
in
t
Healpix_Tables
::
swap_clen
[]
=
const
size_
t
Healpix_Tables
::
swap_clen
[]
=
{
0
,
7
,
5
,
4
,
12
,
10
,
13
,
18
,
14
,
19
,
18
,
17
,
27
,
21
};
const
in
t
Healpix_Tables
::
swap_cycle
[]
=
const
size_
t
Healpix_Tables
::
swap_cycle
[]
=
{
0
,
1
,
8
,
12
,
16
,
21
,
40
,
0
,
1
,
2
,
40
,
114
,
0
,
4
,
160
,
263
,
...
...
Healpix_cxx/healpix_tables.h
View file @
8e9b2689
...
...
@@ -60,7 +60,7 @@ class Healpix_Tables
static
const
int
nb_xoffset
[],
nb_yoffset
[],
nb_facearray
[][
12
],
nb_swaparray
[][
3
];
static
const
in
t
swap_clen
[],
swap_cycle
[];
static
const
size_
t
swap_clen
[],
swap_cycle
[];
};
}}
...
...
mr_util/math_utils.h
View file @
8e9b2689
...
...
@@ -96,7 +96,7 @@ template<typename I> struct isqrt_helper__ <I, true>
static
std
::
uint32_t
isqrt
(
I
arg
)
{
using
namespace
std
;
I
res
=
std
::
sqrt
(
double
(
arg
)
+
0.5
);
I
res
=
I
(
std
::
sqrt
(
double
(
arg
)
+
0.5
)
)
;
if
(
std
::
uint64_t
(
arg
)
<
(
std
::
uint64_t
(
1
)
<<
50
))
return
std
::
uint32_t
(
res
);
if
(
res
*
res
>
arg
)
--
res
;
...
...
@@ -114,7 +114,7 @@ template<typename I> inline std::uint32_t isqrt (I arg)
using
math_utils_detail
::
isqrt
;
/*! Returns the largest integer \a n that fulfills \a 2^n<=arg. */
template
<
typename
I
>
inline
int
ilog2
(
I
arg
)
template
<
typename
I
>
inline
unsigned
int
ilog2
(
I
arg
)
{
#ifdef __GNUC__
if
(
arg
==
0
)
return
0
;
...
...
@@ -125,7 +125,7 @@ template<typename I> inline int ilog2 (I arg)
if
(
sizeof
(
I
)
==
sizeof
(
long
long
))
return
8
*
sizeof
(
long
long
)
-
1
-
__builtin_clzll
(
arg
);
#endif
int
res
=
0
;
unsigned
int
res
=
0
;
while
(
arg
>
0xFFFF
)
{
res
+=
16
;
arg
>>=
16
;
}
if
(
arg
>
0x00FF
)
{
res
|=
8
;
arg
>>=
8
;
}
if
(
arg
>
0x000F
)
{
res
|=
4
;
arg
>>=
4
;
}
...
...
@@ -134,7 +134,7 @@ template<typename I> inline int ilog2 (I arg)
return
res
;
}
template
<
typename
I
>
inline
int
ilog2_nonnull
(
I
arg
)
template
<
typename
I
>
inline
unsigned
int
ilog2_nonnull
(
I
arg
)
{
#ifdef __GNUC__
if
(
sizeof
(
I
)
<=
sizeof
(
int
))
...
...
mr_util/rangeset.h
View file @
8e9b2689
...
...
@@ -88,7 +88,7 @@ template<typename T> class rangeset
size_t
slo
=
sza
<
szb
?
sza
:
szb
,
shi
=
sza
<
szb
?
szb
:
sza
;
double
cost1
=
fct1
*
(
sza
+
szb
);
double
cost2
=
fct2
*
slo
*
std
::
max
(
1
,
ilog2
(
shi
));
double
cost2
=
fct2
*
slo
*
std
::
max
(
1
u
,
ilog2
(
shi
));
return
(
cost1
<=
cost2
)
?
1
:
(
slo
==
sza
)
?
2
:
3
;
}
...
...
Write
Preview
Supports
Markdown
0%
Try again
or
attach a new file
.
Attach a file
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment