Skip to content
GitLab
Projects
Groups
Snippets
Help
Loading...
Help
What's new
10
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Open sidebar
Kai Karius
Fitter
Commits
e8fd0a7f
Commit
e8fd0a7f
authored
Dec 13, 2020
by
karius
Browse files
Options
Browse Files
Download
Email Patches
Plain Diff
Density clean up, added GPUObjects
parent
a85081cc
Changes
13
Hide whitespace changes
Inline
Side-by-side
Showing
13 changed files
with
141 additions
and
682 deletions
+141
-682
Density.cu
Density.cu
+108
-353
Density.h
Density.h
+4
-48
FittingProtocol.h
FittingProtocol.h
+2
-9
Makefile.test
Makefile.test
+3
-0
MemoryResources.h
MemoryResources.h
+12
-7
Parameter.cu
Parameter.cu
+0
-5
Parameter.h
Parameter.h
+11
-1
Query.cu
Query.cu
+0
-78
Query.h
Query.h
+0
-37
Target.cu
Target.cu
+0
-85
Target.h
Target.h
+0
-35
cuda_core.cu
cuda_core.cu
+0
-23
util.h
util.h
+1
-1
No files found.
Density.cu
View file @
e8fd0a7f
...
...
@@ -35,34 +35,6 @@ struct is_outside {
}
};
struct
conditional_to_index
:
public
thrust
::
unary_function
<
thrust
::
tuple
<
uint
,
float
>
,
thrust
::
pair
<
uint3
,
uint3
>>
{
uint
*
d_pixel_dim
;
float
threshold
;
__device__
thrust
::
pair
<
uint3
,
uint3
>
operator
()(
const
thrust
::
tuple
<
uint
,
float
>&
u
)
const
{
thrust
::
pair
<
uint3
,
uint3
>
ret
;
if
(
thrust
::
get
<
1
>
(
u
)
<
threshold
){
//default for min
ret
.
first
.
x
=
UINT32_MAX
;
ret
.
first
.
y
=
UINT32_MAX
;
ret
.
first
.
z
=
UINT32_MAX
;
//default for max
ret
.
second
.
x
=
0
;
ret
.
second
.
y
=
0
;
ret
.
second
.
z
=
0
;
}
else
{
ret
.
first
.
z
=
thrust
::
get
<
0
>
(
u
)
/
(
d_pixel_dim
[
0
]
*
d_pixel_dim
[
1
]);
ret
.
first
.
y
=
(
thrust
::
get
<
0
>
(
u
)
-
ret
.
first
.
z
*
d_pixel_dim
[
0
]
*
d_pixel_dim
[
1
])
/
d_pixel_dim
[
0
];
ret
.
first
.
x
=
thrust
::
get
<
0
>
(
u
)
-
ret
.
first
.
y
*
d_pixel_dim
[
0
]
-
ret
.
first
.
z
*
d_pixel_dim
[
0
]
*
d_pixel_dim
[
1
];
ret
.
second
.
x
=
ret
.
first
.
x
;
ret
.
second
.
y
=
ret
.
first
.
y
;
ret
.
second
.
z
=
ret
.
first
.
z
;
}
return
ret
;
}
};
struct
binarize
:
public
thrust
::
unary_function
<
float
,
float
>
{
float
threshold
;
...
...
@@ -79,17 +51,6 @@ struct binarize : public thrust::unary_function<float,float>
};
struct
uint_bounding_box_pair
:
public
thrust
::
binary_function
<
thrust
::
pair
<
uint3
,
uint3
>
,
thrust
::
pair
<
uint3
,
uint3
>
,
thrust
::
pair
<
uint3
,
uint3
>>
{
__device__
thrust
::
pair
<
uint3
,
uint3
>
operator
()(
const
thrust
::
pair
<
uint3
,
uint3
>&
u
,
const
thrust
::
pair
<
uint3
,
uint3
>&
v
)
const
{
thrust
::
pair
<
uint3
,
uint3
>
ret
;
ret
.
first
=
min
(
u
.
first
,
v
.
first
);
ret
.
second
=
max
(
u
.
second
,
v
.
second
);
return
ret
;
}
};
__global__
void
coords_to_offsets
(
float4
*
d_data
,
uint
*
d_offsets
,
float
*
d_coord_dim
,
uint
*
d_pixel_dim
,
uint
*
d_padding_dim
,
uint
padding_linear_offset
,
size_t
particle_num
){
int
tid
=
blockDim
.
x
*
blockIdx
.
x
+
threadIdx
.
x
;
...
...
@@ -168,45 +129,11 @@ void coords_to_offsets(float4 * d_data, uint * d_offsets, float * d_coord_dim, u
// }
//}
__device__
uint4
Density
::
calc_pixel_indices
(
size_t
pixel_index
){
//TODO: place pixel coords in the middle of the pixel instead of the lower left?
uint4
ret
;
ret
.
z
=
pixel_index
/
(
d_pixel_dim
[
0
]
*
d_pixel_dim
[
1
]);
ret
.
y
=
(
pixel_index
-
ret
.
z
*
d_pixel_dim
[
0
]
*
d_pixel_dim
[
1
])
/
d_pixel_dim
[
1
];
ret
.
x
=
pixel_index
-
ret
.
z
*
d_pixel_dim
[
0
]
*
d_pixel_dim
[
1
]
-
ret
.
y
*
d_pixel_dim
[
1
];
return
ret
;
}
__device__
float
distance
(
const
float4
&
one
,
const
float4
&
other
){
return
sqrtf
(
powf
(
one
.
x
-
other
.
x
,
2.0
)
+
powf
(
one
.
y
-
other
.
y
,
2.0
)
+
powf
(
one
.
z
-
other
.
z
,
2.0
));
}
__device__
bool
Density
::
contains_coords
(
float4
&
coords
){
bool
ret
(
true
);
ret
&=
coords
.
x
<=
d_coord_dim
[
0
];
ret
&=
coords
.
x
>=
0
;
ret
&=
coords
.
y
<=
d_coord_dim
[
1
];
ret
&=
coords
.
y
>=
0
;
ret
&=
coords
.
z
<=
d_coord_dim
[
2
];
ret
&=
coords
.
z
>=
0
;
return
ret
;
}
__device__
bool
Density
::
contains_pixel
(
uint4
&
index
){
bool
ret
(
true
);
ret
&=
d_lower_left
[
0
]
<=
index
.
x
;
ret
&=
d_upper_right
[
0
]
>=
index
.
x
;
ret
&=
d_lower_left
[
1
]
<=
index
.
y
;
ret
&=
d_upper_right
[
1
]
>=
index
.
y
;
ret
&=
d_lower_left
[
2
]
<=
index
.
z
;
ret
&=
d_upper_right
[
2
]
>=
index
.
z
;
return
ret
;
}
__global__
void
min_reduce_strided
(
float
*
d_tmp_array
,
size_t
partition_size
,
int
iteration
){
extern
__shared__
float
mrsshared
[];
...
...
@@ -310,25 +237,15 @@ void copy_result(float * d_tmp_mem, uint * d_indices, size_t index_vol,
distance_map
.
d_data
[
index
]
=
d_tmp_mem
[
tmp_id
];
}
__host__
uint4
Density
::
h_pixel_index_to_pixel
(
const
uint
&
pixel_index
){
uint4
ret
;
ret
.
z
=
pixel_index
/
h_layer_vol
();
ret
.
y
=
(
pixel_index
-
ret
.
z
*
h_layer_vol
())
/
h_bar_vol
();
ret
.
x
=
pixel_index
-
ret
.
z
*
h_layer_vol
()
-
ret
.
y
*
h_bar_vol
();
ret
.
w
=
0
;
return
ret
;
}
__host__
int
*
Density
::
d_above_threshold
(
thrust
::
device_vector
<
int
>
&
above
,
float
t
){
thrust
::
device_ptr
<
float
>
data_ptr
(
d_data
);
is_more_than
pred
{
t
};
thrust
::
fill
(
above
.
begin
(),
above
.
end
(),
0
);
thrust
::
replace_if
(
thrust
::
device
,
above
.
begin
(),
above
.
end
(),
data_ptr
,
pred
,
1
);
int
*
ret
=
thrust
::
raw_pointer_cast
(
&
above
[
0
]);
return
ret
;
}
//__host__
//int * Density::d_above_threshold(thrust::device_vector<int> & above,float t){
// thrust::device_ptr<float> data_ptr(d_data);
// is_more_than pred{t};
// thrust::fill(above.begin(),above.end(),0);
// thrust::replace_if(thrust::device, above.begin(), above.end(), data_ptr, pred, 1);
// int * ret = thrust::raw_pointer_cast(&above[0]);
// return ret;
//}
//__host__
//int * Density::h_above_threshold(thrust::host_vector<int> & above,float t){
...
...
@@ -339,16 +256,6 @@ int * Density::d_above_threshold(thrust::device_vector<int> & above,float t){
// return ret;
//}
__device__
uint4
Density
::
d_pixel_index_to_pixel
(
const
uint
&
pixel_index
){
uint4
ret
;
ret
.
z
=
pixel_index
/
d_layer_vol
();
ret
.
y
=
(
pixel_index
-
ret
.
z
*
d_layer_vol
())
/
d_bar_vol
();
ret
.
x
=
pixel_index
-
ret
.
z
*
d_layer_vol
()
-
ret
.
y
*
d_bar_vol
();
ret
.
w
=
0
;
return
ret
;
}
__device__
float4
Density
::
d_pixel_linear_index_to_coord
(
const
uint
&
pixel_linear_index
){
float4
ret
;
...
...
@@ -490,18 +397,6 @@ struct copy_only_indices : public thrust::unary_function<uint,void>
}
};
__host__
Density
Density
::
only_linear_indeces
(
thrust
::
device_vector
<
uint
>&
td_linear_indices
){
Density
map
;
empty_copy
(
map
,
0.0
);
if
(
td_linear_indices
.
size
()
>
0
){
copy_only_indices
unary_copy_function
;
unary_copy_function
.
d_from
=
d_data
;
unary_copy_function
.
d_to
=
map
.
d_data
;
thrust
::
for_each
(
td_linear_indices
.
begin
(),
td_linear_indices
.
end
(),
unary_copy_function
);
}
return
map
;
}
__global__
void
difference_map_kernel
(
uint
*
d_from
,
uint
*
d_to
,
uint
total_from_chunks
,
uint
to_chunk_size
,
uint
from_chunk_size
,
...
...
@@ -568,13 +463,13 @@ struct linear_offset_to_coord : public thrust::unary_function<uint,float4>
}
};
__host__
std
::
pair
<
float3
,
float3
>
Density
::
get_bounding_box
(
void
){
std
::
pair
<
float3
,
float3
>
ret
;
ret
.
first
=
make_float3
(
h_lower_left
[
0
],
h_lower_left
[
1
],
h_lower_left
[
2
]);
ret
.
second
=
make_float3
(
h_upper_right
[
0
],
h_upper_right
[
1
],
h_upper_right
[
2
]);
return
ret
;
}
//
__host__
//
std::pair<float3,float3> Density::get_bounding_box(void){
//
std::pair<float3,float3> ret;
//
ret.first = make_float3(h_lower_left[0],h_lower_left[1],h_lower_left[2]);
//
ret.second = make_float3(h_upper_right[0],h_upper_right[1],h_upper_right[2]);
//
return ret;
//
}
//__host__
//Particles Density::particles_from_density(const thrust::device_vector<uint>& td_linear_indices){
...
...
@@ -591,29 +486,29 @@ std::pair<float3,float3> Density::get_bounding_box(void){
// return ret;
//}
__host__
void
Density
::
change_canvas
(
uint
*
h_new_pixel_dim
,
uint
*
old_in_new_offset
){
//free and resize device data
cudaMemcpy
(
h_data
,
d_data
,
h_pixel_vol
()
*
sizeof
(
*
h_data
),
cudaMemcpyDeviceToHost
);
cudaFree
(
d_data
);
cudaMalloc
((
void
**
)
&
d_data
,
h_new_pixel_dim
[
0
]
*
h_new_pixel_dim
[
1
]
*
h_new_pixel_dim
[
2
]
*
sizeof
(
*
d_data
));
float
*
h_data_new
=
(
float
*
)
malloc
(
h_new_pixel_dim
[
0
]
*
h_new_pixel_dim
[
1
]
*
h_new_pixel_dim
[
2
]
*
sizeof
(
*
h_data
));
uint3
pixel_index
;
uint
pixel_linear_address
;
for
(
uint
i
=
0
;
i
<
h_pixel_vol
();
i
++
){
pixel_index
=
linear_to_index_space
(
i
,
h_pixel_dim
);
pixel_index
.
x
+=
old_in_new_offset
[
0
];
pixel_index
.
y
+=
old_in_new_offset
[
1
];
pixel_index
.
z
+=
old_in_new_offset
[
2
];
pixel_linear_address
=
index_to_linear_space
(
pixel_index
,
h_new_pixel_dim
);
if
((
pixel_index
.
x
<
h_new_pixel_dim
[
0
])
&&
(
pixel_index
.
y
<
h_new_pixel_dim
[
1
])
&&
(
pixel_index
.
z
<
h_new_pixel_dim
[
2
])){
h_data_new
[
pixel_linear_address
]
=
h_data
[
i
];
}
}
free
(
h_data
);
h_data
=
h_data_new
;
cudaMemcpy
(
d_data
,
h_data
,
h_pixel_vol
()
*
sizeof
(
*
d_data
),
cudaMemcpyHostToDevice
);
}
//
__host__
//
void Density::change_canvas(uint * h_new_pixel_dim, uint * old_in_new_offset){
//
//free and resize device data
//
cudaMemcpy(h_data,d_data,h_pixel_vol()*sizeof(*h_data),cudaMemcpyDeviceToHost);
//
cudaFree(d_data);
//
cudaMalloc((void **) &d_data,h_new_pixel_dim[0]*h_new_pixel_dim[1]*h_new_pixel_dim[2]*sizeof(*d_data));
//
float * h_data_new = (float *) malloc(h_new_pixel_dim[0]*h_new_pixel_dim[1]*h_new_pixel_dim[2]*sizeof(*h_data));
//
uint3 pixel_index;
//
uint pixel_linear_address;
//
for (uint i=0;i<h_pixel_vol();i++){
//
pixel_index = linear_to_index_space(i,h_pixel_dim);
//
pixel_index.x += old_in_new_offset[0];
//
pixel_index.y += old_in_new_offset[1];
//
pixel_index.z += old_in_new_offset[2];
//
pixel_linear_address = index_to_linear_space(pixel_index,h_new_pixel_dim);
//
if ((pixel_index.x < h_new_pixel_dim[0]) && (pixel_index.y < h_new_pixel_dim[1]) && (pixel_index.z < h_new_pixel_dim[2])){
//
h_data_new[pixel_linear_address] = h_data[i];
//
}
//
}
//
free(h_data);
//
h_data = h_data_new;
//
cudaMemcpy(d_data,h_data,h_pixel_vol()*sizeof(*d_data),cudaMemcpyHostToDevice);
//
}
//__host__
//void Density::difference_map_density_threshold(Density & difference_map,float threshold0,float threshold1, thrust::device_vector<uint>& td_on_surface){
...
...
@@ -780,57 +675,6 @@ void Density::to_mrc_async(const char * mrc_path, cudaStream_t &stream){
reader
.
write
();
}
__host__
int
Density
::
add_z_plane
(
int
z
){
int
ret
=
0
;
for
(
int
i
=
0
;
i
<
h_pixel_dim
[
0
];
i
++
){
for
(
int
j
=
0
;
j
<
h_pixel_dim
[
1
];
j
++
){
for
(
int
k
=
0
;
k
<
h_pixel_dim
[
2
];
k
++
){
if
(
k
==
z
)
h_data
[
h_layer_vol
()
*
k
+
h_bar_vol
()
*
j
+
i
]
=
1.0
;
ret
++
;
}
}
}
return
ret
;
}
__host__
int
Density
::
add_x_line
(
int
y
,
int
z
){
int
ret
=
0
;
for
(
int
i
=
0
;
i
<
h_pixel_dim
[
0
];
i
++
){
for
(
int
j
=
0
;
j
<
h_pixel_dim
[
1
];
j
++
){
for
(
int
k
=
0
;
k
<
h_pixel_dim
[
2
];
k
++
){
if
(
j
==
y
&&
k
==
z
)
h_data
[
h_layer_vol
()
*
k
+
h_bar_vol
()
*
j
+
i
]
=
1.0
;
ret
++
;
}
}
}
return
ret
;
}
__host__
int
Density
::
add_relative_sphere
(
const
float
place
[
3
],
float
radius
){
float4
position
;
int
ret
=
0
;
position
.
x
=
h_coord_dim
[
0
]
*
place
[
0
];
position
.
y
=
h_coord_dim
[
1
]
*
place
[
1
];
position
.
z
=
h_coord_dim
[
2
]
*
place
[
2
];
position
.
w
=
0.0
;
radius
*=
h_min_dim
();
float
distance
;
for
(
int
i
=
0
;
i
<
h_pixel_vol
();
i
++
){
distance
=
length
(
h_pixel_linear_index_to_coord
(
i
)
-
position
);
if
(
distance
<=
radius
){
h_data
[
i
]
=
1.0
;
ret
++
;
}
}
return
ret
;
}
//void Density::init_device(){
// if (h_pixel_vol()>0 && host_initiated && !device_initiated){
// if (padding){
...
...
@@ -1132,57 +976,57 @@ void Density::generate_best_n_mrcs(thrust::device_vector<thrust::tuple<float4,fl
nvmlShutdown
();
}
__host__
size_t
Density
::
memory_volume
(){
//in bytes
size_t
vol
=
0
;
//d_coord_dim,3*sizeof(float);
vol
+=
3
*
sizeof
(
float
);
//d_pixel_dim,3*sizeof(uint)
vol
+=
3
*
sizeof
(
uint
);
//d_lower_left,3*sizeof(float)
vol
+=
3
*
sizeof
(
float
);
//d_upper_right,3*sizeof(float)
vol
+=
3
*
sizeof
(
float
);
//d_data
if
(
padding
){
vol
+=
h_padding_dim
[
0
]
*
h_padding_dim
[
1
]
*
h_padding_dim
[
2
]
*
sizeof
(
float
);
}
else
{
vol
+=
h_pixel_vol
()
*
sizeof
(
float
);
}
return
vol
;
}
__host__
void
Density
::
mid_cube_with_surface
(
const
uint3
&
cube_pixel_dim
,
const
uint
&
edge_thickness
,
const
float
&
rho_inner
,
const
float
&
rho_edge
){
uint3
lower_left
=
{
h_pixel_dim
[
0
]
/
2
-
cube_pixel_dim
.
x
/
2
,
h_pixel_dim
[
1
]
/
2
-
cube_pixel_dim
.
y
/
2
,
h_pixel_dim
[
2
]
/
2
-
cube_pixel_dim
.
z
/
2
};
uint3
lower_left_inner
=
{
h_pixel_dim
[
0
]
/
2
-
cube_pixel_dim
.
x
/
2
+
edge_thickness
,
h_pixel_dim
[
1
]
/
2
-
cube_pixel_dim
.
y
/
2
+
edge_thickness
,
h_pixel_dim
[
2
]
/
2
-
cube_pixel_dim
.
z
/
2
+
edge_thickness
};
uint
linear_lower_left
=
index_to_linear_space
(
lower_left
,
h_pixel_dim
);
uint
linear_lower_left_inner
=
index_to_linear_space
(
lower_left_inner
,
h_pixel_dim
);
uint
linear_pixel_offset
;
uint3
index_pixel_offset
;
for
(
uint
z
=
0
;
z
<
cube_pixel_dim
.
z
;
z
++
){
for
(
uint
y
=
0
;
y
<
cube_pixel_dim
.
y
;
y
++
){
for
(
uint
x
=
0
;
x
<
cube_pixel_dim
.
x
;
x
++
){
index_pixel_offset
=
{
x
,
y
,
z
};
linear_pixel_offset
=
linear_lower_left
+
index_to_linear_space
(
index_pixel_offset
,
h_pixel_dim
);
h_data
[
linear_pixel_offset
]
=
rho_edge
;
}
}
}
for
(
uint
z
=
edge_thickness
;
z
<
cube_pixel_dim
.
z
-
edge_thickness
;
z
++
){
for
(
uint
y
=
edge_thickness
;
y
<
cube_pixel_dim
.
y
-
edge_thickness
;
y
++
){
for
(
uint
x
=
edge_thickness
;
x
<
cube_pixel_dim
.
x
-
edge_thickness
;
x
++
){
index_pixel_offset
=
{
x
,
y
,
z
};
linear_pixel_offset
=
linear_lower_left
+
index_to_linear_space
(
index_pixel_offset
,
h_pixel_dim
);
h_data
[
linear_pixel_offset
]
=
rho_inner
;
}
}
}
to_device
();
}
//
__host__
//
size_t Density::memory_volume(){
//
//in bytes
//
size_t vol = 0;
//
//d_coord_dim,3*sizeof(float);
//
vol += 3*sizeof(float);
//
//d_pixel_dim,3*sizeof(uint)
//
vol += 3*sizeof(uint);
//
//d_lower_left,3*sizeof(float)
//
vol += 3*sizeof(float);
//
//d_upper_right,3*sizeof(float)
//
vol += 3*sizeof(float);
//
//d_data
//
if (padding){
//
vol += h_padding_dim[0]*h_padding_dim[1]*h_padding_dim[2]*sizeof(float);
//
} else {
//
vol += h_pixel_vol()*sizeof(float);
//
}
//
return vol;
//
}
//
//
__host__
//
void Density::mid_cube_with_surface(const uint3& cube_pixel_dim, const uint& edge_thickness,
//
const float& rho_inner,const float& rho_edge){
//
uint3 lower_left = {h_pixel_dim[0]/2 - cube_pixel_dim.x/2,h_pixel_dim[1]/2- cube_pixel_dim.y/2,h_pixel_dim[2]/2 - cube_pixel_dim.z/2};
//
uint3 lower_left_inner = {h_pixel_dim[0]/2 - cube_pixel_dim.x/2 + edge_thickness,h_pixel_dim[1]/2- cube_pixel_dim.y/2 + edge_thickness,h_pixel_dim[2]/2 - cube_pixel_dim.z/2 + edge_thickness};
//
uint linear_lower_left = index_to_linear_space(lower_left,h_pixel_dim);
//
uint linear_lower_left_inner = index_to_linear_space(lower_left_inner,h_pixel_dim);
//
uint linear_pixel_offset;
//
uint3 index_pixel_offset;
//
//
for (uint z = 0; z < cube_pixel_dim.z; z++){
//
for (uint y = 0; y < cube_pixel_dim.y; y++){
//
for (uint x=0; x < cube_pixel_dim.x; x++){
//
index_pixel_offset = {x,y,z};
//
linear_pixel_offset = linear_lower_left + index_to_linear_space(index_pixel_offset,h_pixel_dim);
//
h_data[linear_pixel_offset] = rho_edge;
//
}
//
}
//
}
//
for (uint z = edge_thickness; z < cube_pixel_dim.z - edge_thickness; z++){
//
for (uint y = edge_thickness; y < cube_pixel_dim.y - edge_thickness; y++){
//
for (uint x=edge_thickness; x < cube_pixel_dim.x - edge_thickness; x++){
//
index_pixel_offset = {x,y,z};
//
linear_pixel_offset = linear_lower_left +index_to_linear_space(index_pixel_offset,h_pixel_dim);
//
h_data[linear_pixel_offset] = rho_inner;
//
}
//
}
//
}
//
to_device();
//
}
...
...
@@ -1247,50 +1091,9 @@ size_t Density::d_layer_vol(){
return
d_pixel_dim
[
0
]
*
d_pixel_dim
[
1
];
}
__host__
size_t
Density
::
h_bar_vol
(){
return
h_pixel_dim
[
0
];
}
__device__
size_t
Density
::
d_bar_vol
(){
return
d_pixel_dim
[
0
];
}
__host__
float
Density
::
h_min_dim
(){
return
std
::
min
(
h_coord_dim
[
0
],
std
::
min
(
h_coord_dim
[
1
],
h_coord_dim
[
2
]));
}
__device__
float
Density
::
d_min_dim
(){
return
min
(
d_coord_dim
[
0
],
min
(
d_coord_dim
[
1
],
d_coord_dim
[
2
]));
}
__host__
float
Density
::
h_max_dim
(){
return
std
::
max
(
h_coord_dim
[
0
],
std
::
max
(
h_coord_dim
[
1
],
h_coord_dim
[
2
]));
}
__device__
float
Density
::
d_max_dim
(){
return
max
(
d_coord_dim
[
0
],
max
(
d_coord_dim
[
1
],
d_coord_dim
[
2
]));
}
size_t
Density
::
pixel_from_coord_vol
(
float
*
h_label_dim
){
return
vol
(
pixel_from_coord_dim
(
h_label_dim
));
}
__device__
size_t
Density
::
d_linearize_pixel
(
uint3
&
pixel_index
){
return
pixel_index
.
z
*
d_pixel_dim
[
0
]
*
d_pixel_dim
[
1
]
+
pixel_index
.
y
*
d_pixel_dim
[
1
]
+
pixel_index
.
x
;
}
__host__
uint3
Density
::
linear_to_index
_space
(
const
uint
&
linear_index
,
uint
*
const
&
pixel_dim
){
__host__
__device__
uint3
Density
::
linear_to_index
(
const
uint
&
linear_index
,
const
uint3
&
pixel_dim
){
uint3
ret
;
ret
.
z
=
linear_index
/
(
pixel_dim
[
0
]
*
pixel_dim
[
1
]);
ret
.
y
=
(
linear_index
-
ret
.
z
*
pixel_dim
[
0
]
*
pixel_dim
[
1
])
/
pixel_dim
[
0
];
...
...
@@ -1298,6 +1101,20 @@ uint3 Density::linear_to_index_space(const uint &linear_index, uint * const &pix
return
ret
;
}
__host__
__device__
uint
Density
::
index_to_linear
(
const
uint3
&
pixel_index
,
const
uint3
&
pixel_dim
){
return
pixel_index
.
x
+
pixel_index
.
y
*
pixel_dim
.
x
+
pixel_index
.
z
*
pixel_dim
.
y
*
pixel_dim
.
z
;
}
//__host__
//uint3 Density::linear_to_index(const uint &linear_index, uint * const &pixel_dim){
// uint3 ret;
// ret.z = linear_index/(pixel_dim[0]*pixel_dim[1]);
// ret.y = (linear_index - ret.z*pixel_dim[0]*pixel_dim[1])/pixel_dim[0];
// ret.x = linear_index - ret.y*pixel_dim[0] - ret.z*pixel_dim[0]*pixel_dim[1];
// return ret;
//}
__host__
int
Density
::
add_cube_at_with_vol
(
const
uint3
&
pos
,
const
uint3
&
vol
){
int
ret
=
0
;
...
...
@@ -1314,10 +1131,7 @@ int Density::add_cube_at_with_vol(const uint3 &pos, const uint3 &vol){
return
ret
;
}
__host__
uint
Density
::
index_to_linear_space
(
const
uint3
&
pixel_index
,
uint
*
const
&
pixel_dim
){
return
pixel_index
.
x
+
pixel_index
.
y
*
pixel_dim
[
0
]
+
pixel_index
.
z
*
pixel_dim
[
0
]
*
pixel_dim
[
1
];
}
struct
min_functor
:
public
thrust
::
binary_function
<
thrust
::
tuple
<
uint
,
float
>
,
thrust
::
tuple
<
uint
,
float
>
,
thrust
::
tuple
<
uint
,
float
>>
{
...
...
@@ -1485,66 +1299,7 @@ void create_centric_stamp(float *& d_stamp, int *& d_stamp_offsets, const float
// printf("Done\n");
//}
//consumes h_pixel_vol() space
//__host__
//Density Density::density_limits_above(float threshold,uint3 tolerance){
// // typedef these iterators for shorthand
// typedef thrust::device_vector<uint>::iterator UIntIterator;
// typedef thrust::device_vector<float>::iterator FloatIterator;
// // typedef a tuple of these iterators
// typedef thrust::tuple<UIntIterator, FloatIterator> IteratorTuple;
// // typedef the zip_iterator of this tuple
// typedef thrust::zip_iterator<IteratorTuple> ZipIterator;
// // finally, create the zip_iterator
// thrust::device_vector<uint> linear_mem_offsets(h_pixel_vol());
// thrust::sequence(linear_mem_offsets.begin(),linear_mem_offsets.end());
// ZipIterator iter_begin(thrust::make_tuple(linear_mem_offsets.begin(), td_data.begin()));
// ZipIterator iter_end(thrust::make_tuple(linear_mem_offsets.end(), td_data.end()));
// thrust::pair<uint3,uint3> init;
// //default for min
// init.first.x = UINT32_MAX;
// init.first.y = UINT32_MAX;
// init.first.z = UINT32_MAX;
// //default for max
// init.second.x = 0;
// init.second.y = 0;
// init.second.z = 0;
// conditional_to_index to_index;
// to_index.threshold = threshold;
// to_index.d_pixel_dim = d_pixel_dim;
// thrust::pair<uint3,uint3> bounding_box = thrust::transform_reduce(iter_begin, iter_end, to_index, init, uint_bounding_box_pair());
// uint3 index_offset = bounding_box.first-tolerance;
// uint3 pixel_dim = bounding_box.second - bounding_box.first + 2*tolerance;
// uint * h_sub_pixel_dim = (uint *) malloc(3*sizeof(*h_sub_pixel_dim));
// float * h_sub_coord_dim = (float *) malloc(3*sizeof(*h_sub_coord_dim));
// h_sub_pixel_dim[0] = pixel_dim.x;
// h_sub_pixel_dim[1] = pixel_dim.y;
// h_sub_pixel_dim[2] = pixel_dim.z;
//
// for (int i=0;i<3;i++)
// h_sub_coord_dim[i] = (((float) h_sub_pixel_dim[i])/((float) h_pixel_dim[i]))*h_coord_dim[i];
//
// if (pixel_dim.x*pixel_dim.y*pixel_dim.z>0){
// Density sub_density;
// Density::from_bounds(sub_density,h_sub_coord_dim,h_sub_pixel_dim);
// sub_density.fill_with(0.0,true);
// uint3 pixel_index_sub;
// uint3 pixel_index;
// for (uint x = 0; x < h_sub_pixel_dim[0]; x++){
// for (uint y = 0; y < h_sub_pixel_dim[1]; y++){
// for (uint z = 0; z < h_sub_pixel_dim[2]; z++){
// pixel_index = {x + index_offset.x,y + index_offset.y,z + index_offset.z};
// pixel_index_sub = {x,y,z};
// cudaMemcpyAsync(sub_density.d_data + index_to_linear_space(pixel_index_sub,h_sub_pixel_dim),d_data + index_to_linear_space(pixel_index,h_pixel_dim),sizeof(float),cudaMemcpyDeviceToDevice);
// }
// }
// }
// return sub_density;
// } else {
// throw std::length_error("Resulting density seems invalid");
// }
//}
/
//__host__
//void Density::indices_in_range(thrust::device_vector<uint> & td_indices,bool inside, const float & rho0, const float & rho1, size_t & size){
// if (!device_initiated){
...
...
Density.h
View file @
e8fd0a7f
...
...
@@ -35,13 +35,14 @@
#include <chrono>
#include <nvml.h>
#include <GPUObject.h>
#include <MemoryResources.h>
inline
size_t
vol
(
uint4
vec
){
return
vec
.
x
*
vec
.
y
*
vec
.
z
;
}
class
Density
{
class
Density
:
public
GPUObject
{
public:
Density
();
Density
(
float
*
const
&
h_coord_dim
,
const
float
&
h_pixel_size
);
...
...
@@ -70,68 +71,23 @@ public:
//offset meaning: address of first element
uint
padding_linear_offset
{
0
};
uint3
padding_index_offset
=
{
0
,
0
,
0
};
bool
host_initiated
=
false
;
bool
device_initiated
=
false
;
bool
device_updated
=
false
;
void
to_mrc
(
const
char
*
mrc_path
);
void
to_mrc_async
(
const
char
*
mrc_path
,
cudaStream_t
&
stream
);
__host__
size_t
memory_volume
();
__host__
void
change_canvas
(
uint
*
h_new_pixel_dim
,
uint
*
old_in_new_offset
);
__host__
void
mid_cube_with_surface
(
const
uint3
&
cube_pixel_dim
,
const
uint
&
edge_thickness
,
const
float
&
rho_inner
,
const
float
&
rho_edge
);
__host__
std
::
pair
<
float3
,
float3
>
get_bounding_box
(
void
);
__host__
static
void
from_particles_on_grid
(
Density
&
density
,
Particles
&
particles
,
const
float
&
resolution
);
__host__
Density
density_limits_above
(
float
thresold
,
uint3
tolerance
=
{
2
,
2
,
2
});
__host__
void
difference_map_density_threshold
(
Density
&
difference_map
,
float
threshold0
,
float
threshold1
,
thrust
::
device_vector
<
uint
>&
td_on_surface
);
__host__
Density
only_linear_indeces
(
thrust
::
device_vector
<
uint
>&
td_on_surface
);
__host__
void
to_mrc_write_only
(
thrust
::
device_vector
<
uint
>
&
td_linear_indices
,
const
uint
&
label
,
const
char
*
mrc_path
,
uint
log_2_dim
);
__host__
void
cut_and_binarize
(
Density
&
sub_density
,
float
threshold
,
float
above_value
,
float
below_value
,
uint3
tolerance
);
// __host__ Particles particles_from_density(const thrust::device_vector<uint>& td_linear_indices);
__host__
uint3
linear_to_index_space
(
const
uint
&
linear_index
,
uint
*
const
&
pixel_dim
);
__host__
uint
index_to_linear_space
(
const
uint3
&
pixel_index
,
uint
*
const
&
pixel_dim
);
//__host__ uint linear_to_linear_space(const uint &linear_pixel_index, uint * const &pixel_dimA, uint * const &pixel_dimB,uint3 relativ_offset);
__host__
void
generate_best_n_mrcs
(
thrust
::
device_vector
<
thrust
::
tuple
<
float4
,
float3
,
float
>>
td_transformations
);
__host__
int
add_relative_sphere
(
const
float
place
[
3
],
float
radius
);
__host__
int
add_cube_at_with_vol
(
const
uint3
&
pos
,
const
uint3
&
vol
);
__host__
int
add_x_line
(
int
y
,
int
z
);
__host__
int
add_z_plane
(
int
z
);
__host__
int
*
d_above_threshold
(
thrust
::
device_vector
<
int
>
&
above
,
float
t
);
__host__
int
*
h_above_threshold
(
thrust
::
host_vector
<
int
>
&
above
,
float
t
);
__host__
void
fill_with
(
TDensity
,
bool
copy_to_device
=
true
);
__host__
void
empty_copy
(
Density
&
density
,
float
fill
);
__host__
void
init_device
(
void
);
__host__
void
to_device
(
void
);
__host__
void
from_device
();
__host__
void
from_device_async
(
cudaStream_t
&
stream
);
__host__
__device__
static
uint3
linear_to_index
(
const
uint
&
linear_index
,
const
uint3
&
pixel_dim
);
__host__
__device__
static
uint
index_to_linear
(
const
uint3
&
pixel_index
,
const
uint3
&
pixel_dim
);
__host__
size_t
h_pixel_vol
(
void
);
__device__
size_t
d_pixel_vol
(
void
);