# -*- python -*-
#
# Copyright INRIA - CIRAD - INRA
#
# Distributed under the Cecill-C License.
# See accompanying file LICENSE.txt or copy at
# http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html
#
# ==============================================================================
from __future__ import division, print_function
import math
import collections
import cv2
import scipy.spatial
import numpy
import sklearn.neighbors
import openalea.phenomenal.multi_view_reconstruction._c_mvr as c_mvr
from ..object import VoxelGrid
# ==============================================================================
# Class
Voxels = collections.namedtuple("Voxels", ["position", "size"])
VoxelsStage = collections.namedtuple("VoxelsStage", ["consistent", "inconsistent"])
# ==============================================================================
[docs]
def get_voxels_corners(voxels_position, voxels_size):
"""According to the voxels position and their size, return a numpy array
containing for each input voxels the position of the 8 corners.
Parameters
----------
voxels_position : numpy.ndarray
Center position of the voxels
voxels_size : float
Diameter size of the voxels
Returns
-------
a : numpy.array
"""
r = voxels_size / 2.0
x_minus = voxels_position[:, 0] - r
x_plus = voxels_position[:, 0] + r
y_minus = voxels_position[:, 1] - r
y_plus = voxels_position[:, 1] + r
z_minus = voxels_position[:, 2] - r
z_plus = voxels_position[:, 2] + r
a1 = numpy.column_stack((x_minus, y_minus, z_minus))
a2 = numpy.column_stack((x_plus, y_minus, z_minus))
a3 = numpy.column_stack((x_minus, y_plus, z_minus))
a4 = numpy.column_stack((x_minus, y_minus, z_plus))
a5 = numpy.column_stack((x_plus, y_plus, z_minus))
a6 = numpy.column_stack((x_plus, y_minus, z_plus))
a7 = numpy.column_stack((x_minus, y_plus, z_plus))
a8 = numpy.column_stack((x_plus, y_plus, z_plus))
a = numpy.concatenate((a1, a2, a3, a4, a5, a6, a7, a8), axis=1)
a = numpy.reshape(a, (a.shape[0] * 8, 3))
return a
[docs]
def get_bounding_box_voxel_projected(voxels_position, voxels_size, projection):
"""Compute the bounding box value according the radius, angle and
calibration parameters of point_3d projection
Parameters
----------
voxels_position : numpy.ndarray
Center position of voxel
voxels_size : float
Size of side geometry of voxel
projection : function ((x, y, z)) -> (x, y)
Function of projection who take 1 argument (tuple of position (x, y, z))
and return this position 2D (x, y)
Returns
-------
bbox : numpy.ndarray
[[x_min, x_max, y_min, y_max], ...]
Containing min and max value of point_3d projection in x and y axes.
"""
voxels_corners = get_voxels_corners(voxels_position, voxels_size)
pt = projection(voxels_corners)
pt = numpy.reshape(pt, (pt.shape[0] // 8, 8, 2))
bbox = numpy.column_stack((pt.min(axis=1), pt.max(axis=1)))
return bbox
# ==============================================================================
[docs]
def split_voxels_in_eight(voxels):
"""Split each voxel in 8 en return the numpy.array position
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
/ /| / / /|
/ / | /--------/---------/ |
/_ _ _ _ _ _ _ _ _ / | /_ _ _ _ / _ _ _ _ /| |
| | | | | | |/|
| | | =======> | | | / |
| | | |_ _ _ _ | _ _ _ _ |/| |
| | / | | | |/
| | / | | | /
| _ _ _ _ _ _ _ _ _|/ |_ _ _ _ | _ _ _ _ |/
Parameters
----------
voxels : Voxels
Where position is a numpy.array([[x, y, z], ...] containing the center
position of each voxel and size the diameter size value (float) of each
voxel.
Returns
-------
voxels : Voxels
"""
if len(voxels.position) == 0:
return Voxels(voxels.position, voxels.size / 2.0)
r = voxels.size / 4.0
x_minus = voxels.position[:, 0] - r
x_plus = voxels.position[:, 0] + r
y_minus = voxels.position[:, 1] - r
y_plus = voxels.position[:, 1] + r
z_minus = voxels.position[:, 2] - r
z_plus = voxels.position[:, 2] + r
a1 = numpy.column_stack((x_minus, y_minus, z_minus))
a2 = numpy.column_stack((x_plus, y_minus, z_minus))
a3 = numpy.column_stack((x_minus, y_plus, z_minus))
a4 = numpy.column_stack((x_minus, y_minus, z_plus))
a5 = numpy.column_stack((x_plus, y_plus, z_minus))
a6 = numpy.column_stack((x_plus, y_minus, z_plus))
a7 = numpy.column_stack((x_minus, y_plus, z_plus))
a8 = numpy.column_stack((x_plus, y_plus, z_plus))
return Voxels(
numpy.concatenate((a1, a2, a3, a4, a5, a6, a7, a8), axis=0), voxels.size / 2.0
)
# ==============================================================================
[docs]
def voxels_is_visible_in_image(
voxels_position, voxels_size, image, projection, inclusive, image_int=None
):
"""
Return a numpy array containing True if the voxel are
projected is photo-consistent on image else False
**Algorithm**
1. Project each voxel center position on the image, if the position
projected (x, y) is positive on image return True
|
2. If the bounding box of the voxel projected have positive value on
the image the voxel are True
|
3. Check if one pixel containing in the bounding box projected on image
have positive value, if yes return True else return False
Parameters
----------
voxels_position : numpy.array([[x, y, z], ...]
Center position of the voxels
voxels_size : float
diameter size of the voxels
image: numpy.array
Binary image where the voxels are projected.
projection : function (numpy.array([[x, y, z], ...]) -> numpy.array([[x, y], ...])
Function of projection who take 1 argument (numpy.array([[x, y, z],
...] of voxels positions) and return the projected 2D position
numpy.array([[x, y], ...])
inclusive: Describe if the voxels projection are out of the image,
they are considered like still visible
image_int: Integral image of the binary image (optimization)
Returns
-------
out : numpy.array([True, False, ...])
Numpy array containing True if the voxel are
projected is photo-consistent on image else False
"""
height, length = image.shape
ori_result = numpy.zeros(
(
len(
voxels_position,
)
),
dtype=int,
)
r = projection(voxels_position)
cond = (r[:, 0] >= 0) & (r[:, 1] >= 0) & (r[:, 0] < length) & (r[:, 1] < height)
rr = r[cond].astype(int)
(ori_result[cond])[image[rr[:, 1], rr[:, 0]] > 0] = 1
not_cond = numpy.logical_not(ori_result > 0)
result = ori_result[not_cond]
voxels_position = voxels_position[not_cond]
# ==========================================================================
min_xy_max_xy = get_bounding_box_voxel_projected(
voxels_position, voxels_size, projection
)
vv = (
(min_xy_max_xy[:, 2] < 0)
| (min_xy_max_xy[:, 0] >= length)
| (min_xy_max_xy[:, 3] < 0)
| (min_xy_max_xy[:, 1] >= height)
)
not_vv = numpy.logical_not(vv)
result[vv] = 1 if inclusive else 0
min_xy_max_xy = min_xy_max_xy[not_vv]
bb = result[not_vv]
min_xy_max_xy = numpy.floor(min_xy_max_xy).astype(int)
# X limit
(min_xy_max_xy[:, 0])[min_xy_max_xy[:, 0] >= length] = length - 1
(min_xy_max_xy[:, 2])[min_xy_max_xy[:, 2] >= length] = length - 1
# Y limit
(min_xy_max_xy[:, 1])[min_xy_max_xy[:, 1] >= height] = height - 1
(min_xy_max_xy[:, 3])[min_xy_max_xy[:, 3] >= height] = height - 1
# Under zero limit
min_xy_max_xy[:, 0:2] -= 1 # For integral image optimization
min_xy_max_xy[min_xy_max_xy < 0] = 0
# ==========================================================================
for i, (x_min, y_min, x_max, y_max) in enumerate(min_xy_max_xy):
if (
image_int[y_max, x_max]
+ image_int[y_min, x_min]
- image_int[y_min, x_max]
- image_int[y_max, x_min]
) > 0:
bb[i] = 1
# for i, (x_min, y_min, x_max, y_max) in enumerate(min_xy_max_xy):
# if numpy.count_nonzero(image[y_min:y_max + 1, x_min:x_max + 1]) > 0:
# bb[i] = 1
result[not_vv] = bb
ori_result[not_cond] = result
return ori_result
# ==============================================================================
[docs]
def kept_visible_voxel(
voxels_position, voxels_size, image_views, error_tolerance=0, int_images=None
):
"""
Kept in a new collections.deque the voxel who is visible on each image of
images_projections according the error_tolerance
Parameters
----------
voxels_position : numpy.array([[x, y, z], ...]
Center position of the voxels
voxels_size : float
Diameter size of the voxels
image_views : [(image, projection), ...]
List of tuple (image, projection) where image is a binary image
(numpy.ndarray) and function projection (function (x, y, z) -> (x, y))
who take (x, y, z) position on return (x, y) position according space
representation of this image
error_tolerance : int, optional
Number of image will be ignored if the projected voxel is not visible.
int_images: Integral image of the binary image (optimization)
Returns
-------
out : VoxelsStage
"""
photo_consistent = numpy.zeros((len(voxels_position),), dtype=int)
no_kept = None
for i, image_view in enumerate(image_views):
photo_consistent += voxels_is_visible_in_image(
voxels_position,
voxels_size,
image_view.image,
image_view.projection,
image_view.inclusive,
image_int=int_images[i],
)
cond = photo_consistent >= i + 1 - error_tolerance
if no_kept is None:
no_kept = voxels_position[numpy.logical_not(cond)]
else:
no_kept = numpy.insert(
no_kept, 0, voxels_position[numpy.logical_not(cond)], axis=0
)
voxels_position = voxels_position[cond]
photo_consistent = photo_consistent[cond]
consistent = Voxels(voxels_position, voxels_size)
inconsistent = Voxels(no_kept, voxels_size)
return VoxelsStage(consistent, inconsistent)
# ==============================================================================
[docs]
def have_image_ref(image_views):
"""
Returns whether an array of ImageView has a reference image or not.
Parameters
----------
image_views: array[ImageView]
The array of ImageView to test
Returns
-------
True if the array has a reference image else False.
"""
for iv in image_views:
if iv.image_ref is not None:
return True
return False
[docs]
def create_groups(image_views, inconsistent):
groups = collections.defaultdict(list)
kept_groups = collections.defaultdict(list)
group_id = 0
for iv in image_views:
if iv.image_ref is not None:
height, length = iv.image.shape
min_xy_max_xy = get_bounding_box_voxel_projected(
inconsistent.position, inconsistent.size, iv.projection
)
# add each voxel to a visual cones
for i, (x_min, y_min, x_max, y_max) in enumerate(min_xy_max_xy):
if x_max < 0 or x_min >= length or y_max < 0 or y_min >= height:
continue
x_min = int(min(max(math.floor(x_min), 0), length - 1))
y_min = int(min(max(math.floor(y_min), 0), height - 1))
x_max = int(min(max(math.floor(x_max), 0), length - 1))
y_max = int(min(max(math.floor(y_max), 0), height - 1))
img = iv.image_ref[y_min : y_max + 1, x_min : x_max + 1]
yy, xx = numpy.where(img > 0)
yy += y_min
xx += x_min
for y, x in zip(yy, xx):
groups[(group_id, y, x)].append(i)
for y, x in zip(iv.yy, iv.xx):
if len(groups[(group_id, y, x)]) > 0:
kept_groups[(group_id, y, x)] = groups[(group_id, y, x)]
group_id += 1
return kept_groups
[docs]
def check_groups(neigh, inconsistent, groups, nb_distance):
if len(groups.values()) == 0:
return None
positions = []
for index in groups.values():
index = numpy.array(index)
distance, _ = neigh.kneighbors(inconsistent.position[index])
distance = distance.min(axis=1)
xx = distance.argsort()[:nb_distance]
positions.append(inconsistent.position[index[xx]])
position = numpy.unique(numpy.concatenate(positions, axis=0), axis=0)
return Voxels(position, inconsistent.size)
[docs]
def reconstruction_inconsistent(image_views, stages, attractor=None):
for iv in image_views:
if iv.image_ref is not None:
im = project_voxel_centers_on_image(
stages[-1].consistent.position,
stages[-1].consistent.size,
iv.image.shape,
iv.projection,
)
iv.il = iv.image_ref - im
iv.yy, iv.xx = numpy.where(iv.il > 0)
consistent_neighbors = sklearn.neighbors.NearestNeighbors(
n_neighbors=1, metric="euclidean"
)
if numpy.size(stages[-1].consistent.position) == 0:
consistent_neighbors.fit(numpy.array([[0, 0, 0]]))
else:
if attractor is not None:
consistent_neighbors.fit(attractor)
else:
consistent_neighbors.fit(stages[-1].consistent.position)
consistents = [None] * len(stages)
for i, stage in enumerate(stages):
if stage.inconsistent is None:
continue
inconsistent = stage.inconsistent
if consistents[i - 1] is not None:
voxels = split_voxels_in_eight(consistents[i - 1])
position = numpy.concatenate(
(inconsistent.position, voxels.position), axis=0
)
position = numpy.unique(position, axis=0)
inconsistent = Voxels(position, inconsistent.size)
groups = create_groups(image_views, inconsistent)
nb_distance = max(20 - int((20 / len(stages)) * i), 2)
consistents[i] = check_groups(
consistent_neighbors, inconsistent, groups, nb_distance
)
consistent_stages = [None] * len(stages)
for i, (stage, consistent) in enumerate(zip(stages, consistents)):
consistent_stages[i] = stage.consistent
if consistent is not None:
voxels_position = numpy.concatenate(
(consistent_stages[i].position, consistent.position), axis=0
)
voxels_position = numpy.unique(voxels_position, axis=0)
consistent_stages[i] = Voxels(voxels_position, consistent.size)
return consistent_stages
# ==============================================================================
[docs]
def get_integrale_image(img):
a = numpy.zeros_like(img, dtype=int)
a[img > 0] = 1
for y, x in numpy.ndindex(a.shape):
if x - 1 >= 0:
a[y, x] += a[y, x - 1]
if y - 1 >= 0:
a[y, x] += a[y - 1, x]
if x - 1 >= 0 and y - 1 >= 0:
a[y, x] -= a[y - 1, x - 1]
return a
# ==============================================================================
[docs]
def reconstruction_3d(
image_views,
voxels_size=4,
error_tolerance=0,
voxel_center_origin=(0.0, 0.0, 0.0),
start_voxel_size=4096,
voxels_position=None,
attractor=None,
):
"""
Construct a list of voxel represented object with positive value on binary
image in images of images_projections.
Parameters
----------
image_views : [(image, projection), ...]
List of tuple (image, projection) where image is a binary image
(numpy.ndarray) and function projection (function (x, y, z) -> (x, y))
who take (x, y, z) position on return (x, y) position according space
representation of this image
voxels_size : float, optional
Diameter size of the voxels
error_tolerance : int, optional
voxel_center_origin : (x, y, z), optional
Center position of the first original voxel, who will be split.
start_voxel_size: int, optional
Minimum size that the origin voxel size must include at beginning
voxels_position : numpy.ndarray, optional
List of first original voxel who will be split. If None, a list is
created with the voxel_center_origin value.
Returns
-------
out : VoxelGrid
"""
if len(image_views) == 0:
raise ValueError("Len images view have not length")
if voxels_position is None:
voxels_position = numpy.array([voxel_center_origin])
list_voxels_size = [
voxels_size * 2**i
for i in range(20, -1, -1)
if voxels_size * 2 ** (i - 1) < start_voxel_size
]
# Pre-processing (optimization): Compute integral image for speed
# computation
int_images = []
for i, image_view in enumerate(image_views):
a = numpy.zeros_like(image_view.image, dtype=numpy.uint32)
c_mvr.integral_image(image_view.image, a)
int_images.append(a)
stage = VoxelsStage(Voxels(voxels_position, list_voxels_size[0]), None)
stages = [stage]
while stage.consistent.size != voxels_size:
if len(stage.consistent.position) == 0:
break
voxels = split_voxels_in_eight(stage.consistent)
# print(voxels.size)
if voxels.size < 512:
stage = kept_visible_voxel(
voxels.position,
voxels.size,
image_views,
error_tolerance=error_tolerance,
int_images=int_images,
)
else:
stage = VoxelsStage(voxels, None)
stages.append(stage)
consistent_stages = [stage.consistent for stage in stages]
if have_image_ref(image_views):
consistent_stages = reconstruction_inconsistent(
image_views, stages, attractor=attractor
)
return VoxelGrid(consistent_stages[-1].position, consistent_stages[-1].size)
# ==============================================================================
[docs]
def project_voxel_centers_on_image(
voxels_position, voxels_size, shape_image, projection, value=255, dtype=numpy.uint8
):
"""
Create a image with same shape that shape_image and project each voxel on
image and write positive value (255) on it.
Parameters
----------
voxels_position : numpy.ndarray
Voxels center position [[x, y, z], ...]
voxels_size : float
Diameter size of the voxels
shape_image: 2-tuple
Size height and width of the image target projected
projection : function ((x, y, z)) -> (x, y)
Function of projection who take 1 argument (tuple of position (x, y, z))
and return this position 2D (x, y)
value : int
value between 0 and 255 of positive pixel. By default, 255.
dtype : type
numpy type of the returned image. By default, numpy.uint8.
Returns
-------
out : numpy.ndarray
Binary image
"""
height, length = shape_image
img = numpy.zeros((height, length), dtype=dtype)
min_xy_max_xy = get_bounding_box_voxel_projected(
voxels_position, voxels_size, projection
)
vv = (
(min_xy_max_xy[:, 2] < 0)
| (min_xy_max_xy[:, 0] >= length)
| (min_xy_max_xy[:, 3] < 0)
| (min_xy_max_xy[:, 1] >= height)
)
not_vv = numpy.logical_not(vv)
min_xy_max_xy = min_xy_max_xy[not_vv]
min_xy_max_xy = numpy.floor(min_xy_max_xy)
min_xy_max_xy[min_xy_max_xy < 0] = 0
(min_xy_max_xy[:, 0])[min_xy_max_xy[:, 0] >= length] = length - 1
(min_xy_max_xy[:, 1])[min_xy_max_xy[:, 1] >= height] = height - 1
(min_xy_max_xy[:, 2])[min_xy_max_xy[:, 2] >= length] = length - 1
(min_xy_max_xy[:, 3])[min_xy_max_xy[:, 3] >= height] = height - 1
min_xy_max_xy = min_xy_max_xy.astype(int)
for x_min, y_min, x_max, y_max in min_xy_max_xy:
img[y_min : y_max + 1, x_min : x_max + 1] = value
return img
[docs]
def project_voxels_position_on_image(
voxels_position, voxels_size, shape_image, projection
):
"""
Create an image with same shape that shape_image and project each voxel on
image and write positive value (255) on it.
Parameters
----------
voxels_position : [(x, y, z)]
cList (collections.deque) of center position of voxel
voxels_size : float
Size of side geometry of voxel
shape_image: Tuple
size height and length of the image target projected
projection : function ((x, y, z)) -> (x, y)
Function of projection who take 1 argument (tuple of position (x, y, z))
and return this position 2D (x, y)
Returns
-------
out : numpy.ndarray
Binary image
"""
voxels_position = numpy.array(voxels_position)
height, length = shape_image
img = numpy.zeros((height, length), dtype=numpy.uint8)
voxels_corners = get_voxels_corners(voxels_position, voxels_size)
pt = projection(voxels_corners)
pt = numpy.reshape(pt, (pt.shape[0] // 8, 8, 2))
pt[pt < 0] = 0
(pt[:, :, 0])[pt[:, :, 0] >= length] = length - 1
(pt[:, :, 1])[pt[:, :, 1] >= height] = height - 1
pt = numpy.floor(pt).astype(int)
for points in pt:
hull = scipy.spatial.ConvexHull(points)
cv2.fillConvexPoly(img, points[hull.vertices], 255)
return img
# ==============================================================================
[docs]
def image_error(img_ref, img_src, precision=2):
"""
Return false position and true negative result from the comparison on
two binaries images
"""
img_ref = img_ref.astype(numpy.int32)
nb_ref = max(numpy.count_nonzero(img_ref), 1)
nb_ref2 = max(numpy.count_nonzero(img_ref == 0), 1)
img_src = img_src.astype(numpy.int32)
img = numpy.subtract(img_ref, img_src)
true_negative = numpy.bitwise_and(img_ref == 0, img_src == 0)
true_negative = round(
numpy.count_nonzero(true_negative) * 100.0 / nb_ref2, precision
)
# true_negative = numpy.bitwise_and(img_ref == 0, img_src == 0)
false_positive = round(
numpy.count_nonzero(img[img < 0]) * 100.0 / nb_ref2, precision
)
false_negative = round(
numpy.count_nonzero(img[img > 0]) * 100.0 / nb_ref, precision
)
print(true_negative, false_positive, false_negative)
return false_positive, false_negative
[docs]
def reconstruction_error(voxels_grid, image_views):
"""
Compute the reconstruction error (false positive and true negative) of
the 3d reconstruction from the image view.
Parameters
----------
voxels_grid: VoxelGrid
The voxel grid
image_views : numpy.ndarray[ImageView]
An array of all image views
Returns
-------
out : (float,float)
A tuple with the mean false positive and the mean false negative
"""
sum_false_positive = 0
sum_false_negative = 0
for image_view in image_views:
img_src = project_voxel_centers_on_image(
voxels_grid.voxels_position,
voxels_grid.voxels_size,
image_view.image.shape,
image_view.projection,
)
false_positive, false_negative = image_error(image_view.image, img_src)
sum_false_positive += false_positive
sum_false_negative += false_negative
mean_false_positive = sum_false_positive / len(image_views)
mean_false_negative = sum_false_negative / len(image_views)
return mean_false_positive, mean_false_negative
