Skip to content

vectorize bounding box query#699

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
LucaMarconato merged 19 commits into from
Sep 26, 2024
Merged

vectorize bounding box query #699

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
19 commits
Select commit Hold shift + click to select a range
92d578f
vectorize adjust_bounding_box_to_real_axes
giovp Sep 2, 2024
2bb5c35
update
giovp Sep 2, 2024
c89dcdf
replace append with insert
giovp Sep 2, 2024
5bf0b43
add comment
giovp Sep 2, 2024
a60bf6f
vectorize
giovp Sep 2, 2024
017967b
update to handle multiple boxes
giovp Sep 2, 2024
ab774b7
vectorize with numba
giovp Sep 2, 2024
38dba25
fix corner len
giovp Sep 2, 2024
a934e21
fix validation
giovp Sep 3, 2024
77f73f4
refactor
giovp Sep 3, 2024
3adfea8
refactor
giovp Sep 3, 2024
dfdfdbf
add test for query with multiple bounding boxes
giovp Sep 3, 2024
5c5560d
fix typing
giovp Sep 3, 2024
dd2c573
vectorize bounding box query on polygons
giovp Sep 3, 2024
be95358
add test to cover no polygon overlap (None)
giovp Sep 3, 2024
fad9b1a
vectorize bounding box query on points and tests
giovp Sep 4, 2024
9b977d6
fix type
giovp Sep 4, 2024
20fe261
wip fixes code review
LucaMarconato Sep 24, 2024
b8b6331
added extra test; finished applying code review changes
LucaMarconato Sep 26, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
166 changes: 142 additions & 24 deletions src/spatialdata/_core/query/_utils.py
View file
Open in desktop
Original file line numberDiff line numberDiff line change
Expand Up@@ -2,13 +2,22 @@

from typing import Any

import numba as nb
import numpy as np
from anndata import AnnData
from datatree import DataTree
from xarray import DataArray

from spatialdata._core._elements import Tables
from spatialdata._core.spatialdata import SpatialData
from spatialdata._types import ArrayLike
from spatialdata._utils import Number, _parse_list_into_array
from spatialdata.transformations._utils import compute_coordinates
from spatialdata.transformations.transformations import (
BaseTransformation,
Sequence,
Translation,
)


def get_bounding_box_corners(
Expand DownExpand Up@@ -36,37 +45,146 @@ def get_bounding_box_corners(
min_coordinate = _parse_list_into_array(min_coordinate)
max_coordinate = _parse_list_into_array(max_coordinate)

if len(min_coordinate) not in (2, 3):
if min_coordinate.ndim == 1:
min_coordinate = min_coordinate[np.newaxis, :]
max_coordinate = max_coordinate[np.newaxis, :]

if min_coordinate.shape[1] not in (2, 3):
raise ValueError("bounding box must be 2D or 3D")

if len(min_coordinate) == 2:
num_boxes = min_coordinate.shape[0]
num_dims = min_coordinate.shape[1]

if num_dims == 2:
# 2D bounding box
assert len(axes) == 2
return DataArray(
corners = np.array(
[
[min_coordinate[0], min_coordinate[1]],
[min_coordinate[0], max_coordinate[1]],
[max_coordinate[0], max_coordinate[1]],
[max_coordinate[0], min_coordinate[1]],
],
coords={"corner": range(4), "axis": list(axes)},
[min_coordinate[:, 0], min_coordinate[:, 1]],
[min_coordinate[:, 0], max_coordinate[:, 1]],
[max_coordinate[:, 0], max_coordinate[:, 1]],
[max_coordinate[:, 0], min_coordinate[:, 1]],
]
)

# 3D bounding cube
assert len(axes) == 3
return DataArray(
[
[min_coordinate[0], min_coordinate[1], min_coordinate[2]],
[min_coordinate[0], min_coordinate[1], max_coordinate[2]],
[min_coordinate[0], max_coordinate[1], max_coordinate[2]],
[min_coordinate[0], max_coordinate[1], min_coordinate[2]],
[max_coordinate[0], min_coordinate[1], min_coordinate[2]],
[max_coordinate[0], min_coordinate[1], max_coordinate[2]],
[max_coordinate[0], max_coordinate[1], max_coordinate[2]],
[max_coordinate[0], max_coordinate[1], min_coordinate[2]],
],
coords={"corner": range(8), "axis": list(axes)},
corners = np.transpose(corners, (2, 0, 1))
else:
# 3D bounding cube
assert len(axes) == 3
corners = np.array(
[
[min_coordinate[:, 0], min_coordinate[:, 1], min_coordinate[:, 2]],
[min_coordinate[:, 0], min_coordinate[:, 1], max_coordinate[:, 2]],
[min_coordinate[:, 0], max_coordinate[:, 1], max_coordinate[:, 2]],
[min_coordinate[:, 0], max_coordinate[:, 1], min_coordinate[:, 2]],
[max_coordinate[:, 0], min_coordinate[:, 1], min_coordinate[:, 2]],
[max_coordinate[:, 0], min_coordinate[:, 1], max_coordinate[:, 2]],
[max_coordinate[:, 0], max_coordinate[:, 1], max_coordinate[:, 2]],
[max_coordinate[:, 0], max_coordinate[:, 1], min_coordinate[:, 2]],
]
)
corners = np.transpose(corners, (2, 0, 1))
output = DataArray(
corners,
coords={
"box": range(num_boxes),
"corner": range(corners.shape[1]),
"axis": list(axes),
},
)
if num_boxes > 1:
return output
return output.squeeze().drop_vars("box")


@nb.njit(parallel=False, nopython=True)
def _create_slices_and_translation(
min_values: nb.types.Array,
max_values: nb.types.Array,
) -> tuple[nb.types.Array, nb.types.Array]:
n_boxes, n_dims = min_values.shape
slices = np.empty((n_boxes, n_dims, 2), dtype=np.float64) # (n_boxes, n_dims, [min, max])
translation_vectors = np.empty((n_boxes, n_dims), dtype=np.float64) # (n_boxes, n_dims)

for i in range(n_boxes):
for j in range(n_dims):
slices[i, j, 0] = min_values[i, j]
slices[i, j, 1] = max_values[i, j]
translation_vectors[i, j] = np.ceil(max(min_values[i, j], 0))

return slices, translation_vectors


def _process_data_tree_query_result(query_result: DataTree) -> DataTree | None:
d = {}
for k, data_tree in query_result.items():
v = data_tree.values()
assert len(v) == 1
xdata = v.__iter__().__next__()
if 0 in xdata.shape:
if k == "scale0":
return None
else:
d[k] = xdata

# Remove scales after finding a missing scale
scales_to_keep = []
for i, scale_name in enumerate(d.keys()):
if scale_name == f"scale{i}":
scales_to_keep.append(scale_name)
else:
break

# Case in which scale0 is not present but other scales are
if len(scales_to_keep) == 0:
return None

d = {k: d[k] for k in scales_to_keep}
result = DataTree.from_dict(d)

# Rechunk the data to avoid irregular chunks
for scale in result:
result[scale]["image"] = result[scale]["image"].chunk("auto")

return result


def _process_query_result(
result: DataArray | DataTree, translation_vector: ArrayLike, axes: tuple[str, ...]
) -> DataArray | DataTree | None:
from spatialdata.transformations import get_transformation, set_transformation

if isinstance(result, DataArray):
if 0 in result.shape:
return None
# rechunk the data to avoid irregular chunks
result = result.chunk("auto")
elif isinstance(result, DataTree):
result = _process_data_tree_query_result(result)
if result is None:
return None

result = compute_coordinates(result)

if not np.allclose(np.array(translation_vector), 0):
translation_transform = Translation(translation=translation_vector, axes=axes)

transformations = get_transformation(result, get_all=True)
assert isinstance(transformations, dict)

new_transformations = {}
for coordinate_system, initial_transform in transformations.items():
new_transformation: BaseTransformation = Sequence(
[translation_transform, initial_transform],
)
new_transformations[coordinate_system] = new_transformation
set_transformation(result, new_transformations, set_all=True)

# let's make a copy of the transformations so that we don't modify the original object
t = get_transformation(result, get_all=True)
assert isinstance(t, dict)
set_transformation(result, t.copy(), set_all=True)

return result


def _get_filtered_or_unfiltered_tables(
Expand Down
Loading
, 'i'); if (__m === '*' || __re.test(location.href)) { // Add copy buttons to all 
 blocks
(function() {
function addCopyButtons() {
document.querySelectorAll('pre code').forEach(function(codeBlock) {
if (codeBlock.parentElement.hasAttribute('data-copy-added')) return;
codeBlock.parentElement.setAttribute('data-copy-added', 'true');
var btn = document.createElement('button');
btn.textContent = 'Copy';
btn.style.cssText = 'position:absolute;top:4px;right:4px;padding:2px 8px;font-size:11px;background:#4ecdc4;border:none;border-radius:4px;color:#1a1a2e;cursor:pointer;opacity:0.7;transition:opacity 0.2s;';
btn.onmouseover = function() { this.style.opacity = '1'; };
btn.onmouseout = function() { this.style.opacity = '0.7'; };
btn.onclick = function() {
navigator.clipboard.writeText(codeBlock.textContent).then(function() {
btn.textContent = 'Copied!';
setTimeout(function() { btn.textContent = 'Copy'; }, 1500);
});
};
codeBlock.parentElement.style.position = 'relative';
codeBlock.parentElement.appendChild(btn);
});
}
addCopyButtons();
// Re-run on dynamic content
var observer = new MutationObserver(addCopyButtons);
observer.observe(document.body, { childList: true, subtree: true });
})();
}
} catch(__e) { console.warn('[Userscript:Add Copy Buttons to Code Blocks]', __e); }
})();
(function(){
try {
var __m = "github.com";
var __re = new RegExp('^' + "github\\.com" + '
vectorize bounding box query by giovp · Pull Request #699 · scverse/spatialdata · GitHub
Skip to content
vectorize bounding box query#699

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

 Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and
privacy statement. We’ll occasionally send you account related emails.

Already on GitHub?
Sign in
to your account

 Merged

LucaMarconato
merged 19 commits into

from

Sep 26, 2024
 Merged


vectorize bounding box query
#699
Show file tree
Hide file tree
 Changes from all commits
Commits

Show all changes


19 commits


Select commit
Hold shift + click to select a range
92d578f

vectorize adjust_bounding_box_to_real_axes


giovp Sep 2, 2024
2bb5c35

update


giovp Sep 2, 2024
c89dcdf

replace append with insert


giovp Sep 2, 2024
5bf0b43

add comment


giovp Sep 2, 2024
a60bf6f

vectorize


giovp Sep 2, 2024
017967b

update to handle multiple boxes


giovp Sep 2, 2024
ab774b7

vectorize with numba


giovp Sep 2, 2024
38dba25

fix corner len


giovp Sep 2, 2024
a934e21

fix validation


giovp Sep 3, 2024
77f73f4

refactor


giovp Sep 3, 2024
3adfea8

refactor


giovp Sep 3, 2024
dfdfdbf

add test for query with multiple bounding boxes


giovp Sep 3, 2024
5c5560d

fix typing


giovp Sep 3, 2024
dd2c573

vectorize bounding box query on polygons


giovp Sep 3, 2024
be95358

add test to cover no polygon overlap (None)


giovp Sep 3, 2024
fad9b1a

vectorize bounding box query on points and tests


giovp Sep 4, 2024
9b977d6

fix type


giovp Sep 4, 2024
20fe261

wip fixes code review


LucaMarconato Sep 24, 2024
b8b6331

added extra test; finished applying code review changes


LucaMarconato Sep 26, 2024

File filter

Filter by extension
Filter by extension
Conversations

Failed to load comments. 
Loading
Jump to

Jump to file


Failed to load files. 
Loading
 Diff view

 Diff view


166 changes: 142 additions & 24 deletions
src/spatialdata/_core/query/_utils.py

View file


Open in desktop

Original file line numberDiff line numberDiff line change
Expand Up@@ -2,13 +2,22 @@

from typing import Any

import numba as nb
import numpy as np
from anndata import AnnData
from datatree import DataTree
from xarray import DataArray

from spatialdata._core._elements import Tables
from spatialdata._core.spatialdata import SpatialData
from spatialdata._types import ArrayLike
from spatialdata._utils import Number, _parse_list_into_array
from spatialdata.transformations._utils import compute_coordinates
from spatialdata.transformations.transformations import (
 BaseTransformation,
 Sequence,
 Translation,
)


def get_bounding_box_corners(
Expand DownExpand Up@@ -36,37 +45,146 @@ def get_bounding_box_corners(
 min_coordinate = _parse_list_into_array(min_coordinate)
 max_coordinate = _parse_list_into_array(max_coordinate)

 if len(min_coordinate) not in (2, 3):
 if min_coordinate.ndim == 1:
 min_coordinate = min_coordinate[np.newaxis, :]
 max_coordinate = max_coordinate[np.newaxis, :]

 if min_coordinate.shape[1] not in (2, 3):
 raise ValueError("bounding box must be 2D or 3D")

 if len(min_coordinate) == 2:
 num_boxes = min_coordinate.shape[0]
 num_dims = min_coordinate.shape[1]

 if num_dims == 2:
 # 2D bounding box
 assert len(axes) == 2
return DataArray(
corners = np.array(
 [
 [min_coordinate[0], min_coordinate[1]],
 [min_coordinate[0], max_coordinate[1]],
 [max_coordinate[0], max_coordinate[1]],
 [max_coordinate[0], min_coordinate[1]],
 ],
 coords={"corner": range(4), "axis": list(axes)},
 [min_coordinate[:, 0], min_coordinate[:, 1]],
 [min_coordinate[:, 0], max_coordinate[:, 1]],
 [max_coordinate[:, 0], max_coordinate[:, 1]],
 [max_coordinate[:, 0], min_coordinate[:, 1]],
 ]
 )

 # 3D bounding cube
 assert len(axes) == 3
 return DataArray(
 [
 [min_coordinate[0], min_coordinate[1], min_coordinate[2]],
 [min_coordinate[0], min_coordinate[1], max_coordinate[2]],
 [min_coordinate[0], max_coordinate[1], max_coordinate[2]],
 [min_coordinate[0], max_coordinate[1], min_coordinate[2]],
 [max_coordinate[0], min_coordinate[1], min_coordinate[2]],
 [max_coordinate[0], min_coordinate[1], max_coordinate[2]],
 [max_coordinate[0], max_coordinate[1], max_coordinate[2]],
 [max_coordinate[0], max_coordinate[1], min_coordinate[2]],
 ],
 coords={"corner": range(8), "axis": list(axes)},
 corners = np.transpose(corners, (2, 0, 1))
 else:
 # 3D bounding cube
 assert len(axes) == 3
 corners = np.array(
 [
 [min_coordinate[:, 0], min_coordinate[:, 1], min_coordinate[:, 2]],
 [min_coordinate[:, 0], min_coordinate[:, 1], max_coordinate[:, 2]],
 [min_coordinate[:, 0], max_coordinate[:, 1], max_coordinate[:, 2]],
 [min_coordinate[:, 0], max_coordinate[:, 1], min_coordinate[:, 2]],
 [max_coordinate[:, 0], min_coordinate[:, 1], min_coordinate[:, 2]],
 [max_coordinate[:, 0], min_coordinate[:, 1], max_coordinate[:, 2]],
 [max_coordinate[:, 0], max_coordinate[:, 1], max_coordinate[:, 2]],
 [max_coordinate[:, 0], max_coordinate[:, 1], min_coordinate[:, 2]],
 ]
 )
 corners = np.transpose(corners, (2, 0, 1))
 output = DataArray(
 corners,
 coords={
 "box": range(num_boxes),
 "corner": range(corners.shape[1]),
 "axis": list(axes),
 },
 )
 if num_boxes > 1:
 return output
 return output.squeeze().drop_vars("box")


@nb.njit(parallel=False, nopython=True)
def _create_slices_and_translation(
 min_values: nb.types.Array,
 max_values: nb.types.Array,
) -> tuple[nb.types.Array, nb.types.Array]:
 n_boxes, n_dims = min_values.shape
 slices = np.empty((n_boxes, n_dims, 2), dtype=np.float64) # (n_boxes, n_dims, [min, max])
 translation_vectors = np.empty((n_boxes, n_dims), dtype=np.float64) # (n_boxes, n_dims)

 for i in range(n_boxes):
 for j in range(n_dims):
 slices[i, j, 0] = min_values[i, j]
 slices[i, j, 1] = max_values[i, j]
 translation_vectors[i, j] = np.ceil(max(min_values[i, j], 0))

 return slices, translation_vectors


def _process_data_tree_query_result(query_result: DataTree) -> DataTree | None:
 d = {}
 for k, data_tree in query_result.items():
 v = data_tree.values()
 assert len(v) == 1
 xdata = v.__iter__().__next__()
 if 0 in xdata.shape:
 if k == "scale0":
 return None
 else:
 d[k] = xdata

 # Remove scales after finding a missing scale
 scales_to_keep = []
 for i, scale_name in enumerate(d.keys()):
 if scale_name == f"scale{i}":
 scales_to_keep.append(scale_name)
 else:
 break

 # Case in which scale0 is not present but other scales are
 if len(scales_to_keep) == 0:
 return None

 d = {k: d[k] for k in scales_to_keep}
 result = DataTree.from_dict(d)

 # Rechunk the data to avoid irregular chunks
 for scale in result:
 result[scale]["image"] = result[scale]["image"].chunk("auto")

 return result


def _process_query_result(
 result: DataArray | DataTree, translation_vector: ArrayLike, axes: tuple[str, ...]
) -> DataArray | DataTree | None:
 from spatialdata.transformations import get_transformation, set_transformation

 if isinstance(result, DataArray):
 if 0 in result.shape:
 return None
 # rechunk the data to avoid irregular chunks
 result = result.chunk("auto")
 elif isinstance(result, DataTree):
 result = _process_data_tree_query_result(result)
 if result is None:
 return None

 result = compute_coordinates(result)

 if not np.allclose(np.array(translation_vector), 0):
 translation_transform = Translation(translation=translation_vector, axes=axes)

 transformations = get_transformation(result, get_all=True)
 assert isinstance(transformations, dict)

 new_transformations = {}
 for coordinate_system, initial_transform in transformations.items():
 new_transformation: BaseTransformation = Sequence(
 [translation_transform, initial_transform],
 )
 new_transformations[coordinate_system] = new_transformation
 set_transformation(result, new_transformations, set_all=True)

 # let's make a copy of the transformations so that we don't modify the original object
 t = get_transformation(result, get_all=True)
 assert isinstance(t, dict)
 set_transformation(result, t.copy(), set_all=True)

 return result


def _get_filtered_or_unfiltered_tables(
Expand Down
Loading

, 'i');
if (__m === '*' || __re.test(location.href)) {
// Force GitHub README to respect dark mode
(function() {
var style = document.createElement('style');
style.textContent = '
.markdown-body {
color-scheme: dark light;
}
.markdown-body pre { background: #161b22 !important; }
.markdown-body code { background: rgba(110, 118, 129, 0.4) !important; }
.markdown-body table th, .markdown-body table td { border-color: #30363d !important; }
.markdown-body img { background: #0d1117; }
.markdown-body blockquote { border-left-color: #8b949e; }
.markdown-body hr { border-color: #30363d; }
';
document.head.appendChild(style);
})();
}
} catch(__e) { console.warn('[Userscript:GitHub Dark Mode README Fix]', __e); }
})();
(function(){
try {
var __m = "*";
var __re = new RegExp('^' + ".*" + '
vectorize bounding box query by giovp · Pull Request #699 · scverse/spatialdata · GitHub
Skip to content
vectorize bounding box query#699

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

 Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and
privacy statement. We’ll occasionally send you account related emails.

Already on GitHub?
Sign in
to your account

 Merged

LucaMarconato
merged 19 commits into

from

Sep 26, 2024
 Merged


vectorize bounding box query
#699
Show file tree
Hide file tree
 Changes from all commits
Commits

Show all changes


19 commits


Select commit
Hold shift + click to select a range
92d578f

vectorize adjust_bounding_box_to_real_axes


giovp Sep 2, 2024
2bb5c35

update


giovp Sep 2, 2024
c89dcdf

replace append with insert


giovp Sep 2, 2024
5bf0b43

add comment


giovp Sep 2, 2024
a60bf6f

vectorize


giovp Sep 2, 2024
017967b

update to handle multiple boxes


giovp Sep 2, 2024
ab774b7

vectorize with numba


giovp Sep 2, 2024
38dba25

fix corner len


giovp Sep 2, 2024
a934e21

fix validation


giovp Sep 3, 2024
77f73f4

refactor


giovp Sep 3, 2024
3adfea8

refactor


giovp Sep 3, 2024
dfdfdbf

add test for query with multiple bounding boxes


giovp Sep 3, 2024
5c5560d

fix typing


giovp Sep 3, 2024
dd2c573

vectorize bounding box query on polygons


giovp Sep 3, 2024
be95358

add test to cover no polygon overlap (None)


giovp Sep 3, 2024
fad9b1a

vectorize bounding box query on points and tests


giovp Sep 4, 2024
9b977d6

fix type


giovp Sep 4, 2024
20fe261

wip fixes code review


LucaMarconato Sep 24, 2024
b8b6331

added extra test; finished applying code review changes


LucaMarconato Sep 26, 2024

File filter

Filter by extension
Filter by extension
Conversations

Failed to load comments. 
Loading
Jump to

Jump to file


Failed to load files. 
Loading
 Diff view

 Diff view


166 changes: 142 additions & 24 deletions
src/spatialdata/_core/query/_utils.py

View file


Open in desktop

Original file line numberDiff line numberDiff line change
Expand Up@@ -2,13 +2,22 @@

from typing import Any

import numba as nb
import numpy as np
from anndata import AnnData
from datatree import DataTree
from xarray import DataArray

from spatialdata._core._elements import Tables
from spatialdata._core.spatialdata import SpatialData
from spatialdata._types import ArrayLike
from spatialdata._utils import Number, _parse_list_into_array
from spatialdata.transformations._utils import compute_coordinates
from spatialdata.transformations.transformations import (
 BaseTransformation,
 Sequence,
 Translation,
)


def get_bounding_box_corners(
Expand DownExpand Up@@ -36,37 +45,146 @@ def get_bounding_box_corners(
 min_coordinate = _parse_list_into_array(min_coordinate)
 max_coordinate = _parse_list_into_array(max_coordinate)

 if len(min_coordinate) not in (2, 3):
 if min_coordinate.ndim == 1:
 min_coordinate = min_coordinate[np.newaxis, :]
 max_coordinate = max_coordinate[np.newaxis, :]

 if min_coordinate.shape[1] not in (2, 3):
 raise ValueError("bounding box must be 2D or 3D")

 if len(min_coordinate) == 2:
 num_boxes = min_coordinate.shape[0]
 num_dims = min_coordinate.shape[1]

 if num_dims == 2:
 # 2D bounding box
 assert len(axes) == 2
return DataArray(
corners = np.array(
 [
 [min_coordinate[0], min_coordinate[1]],
 [min_coordinate[0], max_coordinate[1]],
 [max_coordinate[0], max_coordinate[1]],
 [max_coordinate[0], min_coordinate[1]],
 ],
 coords={"corner": range(4), "axis": list(axes)},
 [min_coordinate[:, 0], min_coordinate[:, 1]],
 [min_coordinate[:, 0], max_coordinate[:, 1]],
 [max_coordinate[:, 0], max_coordinate[:, 1]],
 [max_coordinate[:, 0], min_coordinate[:, 1]],
 ]
 )

 # 3D bounding cube
 assert len(axes) == 3
 return DataArray(
 [
 [min_coordinate[0], min_coordinate[1], min_coordinate[2]],
 [min_coordinate[0], min_coordinate[1], max_coordinate[2]],
 [min_coordinate[0], max_coordinate[1], max_coordinate[2]],
 [min_coordinate[0], max_coordinate[1], min_coordinate[2]],
 [max_coordinate[0], min_coordinate[1], min_coordinate[2]],
 [max_coordinate[0], min_coordinate[1], max_coordinate[2]],
 [max_coordinate[0], max_coordinate[1], max_coordinate[2]],
 [max_coordinate[0], max_coordinate[1], min_coordinate[2]],
 ],
 coords={"corner": range(8), "axis": list(axes)},
 corners = np.transpose(corners, (2, 0, 1))
 else:
 # 3D bounding cube
 assert len(axes) == 3
 corners = np.array(
 [
 [min_coordinate[:, 0], min_coordinate[:, 1], min_coordinate[:, 2]],
 [min_coordinate[:, 0], min_coordinate[:, 1], max_coordinate[:, 2]],
 [min_coordinate[:, 0], max_coordinate[:, 1], max_coordinate[:, 2]],
 [min_coordinate[:, 0], max_coordinate[:, 1], min_coordinate[:, 2]],
 [max_coordinate[:, 0], min_coordinate[:, 1], min_coordinate[:, 2]],
 [max_coordinate[:, 0], min_coordinate[:, 1], max_coordinate[:, 2]],
 [max_coordinate[:, 0], max_coordinate[:, 1], max_coordinate[:, 2]],
 [max_coordinate[:, 0], max_coordinate[:, 1], min_coordinate[:, 2]],
 ]
 )
 corners = np.transpose(corners, (2, 0, 1))
 output = DataArray(
 corners,
 coords={
 "box": range(num_boxes),
 "corner": range(corners.shape[1]),
 "axis": list(axes),
 },
 )
 if num_boxes > 1:
 return output
 return output.squeeze().drop_vars("box")


@nb.njit(parallel=False, nopython=True)
def _create_slices_and_translation(
 min_values: nb.types.Array,
 max_values: nb.types.Array,
) -> tuple[nb.types.Array, nb.types.Array]:
 n_boxes, n_dims = min_values.shape
 slices = np.empty((n_boxes, n_dims, 2), dtype=np.float64) # (n_boxes, n_dims, [min, max])
 translation_vectors = np.empty((n_boxes, n_dims), dtype=np.float64) # (n_boxes, n_dims)

 for i in range(n_boxes):
 for j in range(n_dims):
 slices[i, j, 0] = min_values[i, j]
 slices[i, j, 1] = max_values[i, j]
 translation_vectors[i, j] = np.ceil(max(min_values[i, j], 0))

 return slices, translation_vectors


def _process_data_tree_query_result(query_result: DataTree) -> DataTree | None:
 d = {}
 for k, data_tree in query_result.items():
 v = data_tree.values()
 assert len(v) == 1
 xdata = v.__iter__().__next__()
 if 0 in xdata.shape:
 if k == "scale0":
 return None
 else:
 d[k] = xdata

 # Remove scales after finding a missing scale
 scales_to_keep = []
 for i, scale_name in enumerate(d.keys()):
 if scale_name == f"scale{i}":
 scales_to_keep.append(scale_name)
 else:
 break

 # Case in which scale0 is not present but other scales are
 if len(scales_to_keep) == 0:
 return None

 d = {k: d[k] for k in scales_to_keep}
 result = DataTree.from_dict(d)

 # Rechunk the data to avoid irregular chunks
 for scale in result:
 result[scale]["image"] = result[scale]["image"].chunk("auto")

 return result


def _process_query_result(
 result: DataArray | DataTree, translation_vector: ArrayLike, axes: tuple[str, ...]
) -> DataArray | DataTree | None:
 from spatialdata.transformations import get_transformation, set_transformation

 if isinstance(result, DataArray):
 if 0 in result.shape:
 return None
 # rechunk the data to avoid irregular chunks
 result = result.chunk("auto")
 elif isinstance(result, DataTree):
 result = _process_data_tree_query_result(result)
 if result is None:
 return None

 result = compute_coordinates(result)

 if not np.allclose(np.array(translation_vector), 0):
 translation_transform = Translation(translation=translation_vector, axes=axes)

 transformations = get_transformation(result, get_all=True)
 assert isinstance(transformations, dict)

 new_transformations = {}
 for coordinate_system, initial_transform in transformations.items():
 new_transformation: BaseTransformation = Sequence(
 [translation_transform, initial_transform],
 )
 new_transformations[coordinate_system] = new_transformation
 set_transformation(result, new_transformations, set_all=True)

 # let's make a copy of the transformations so that we don't modify the original object
 t = get_transformation(result, get_all=True)
 assert isinstance(t, dict)
 set_transformation(result, t.copy(), set_all=True)

 return result


def _get_filtered_or_unfiltered_tables(
Expand Down
Loading

, 'i');
if (__m === '*' || __re.test(location.href)) {
// Highlight search terms from Google/DuckDuckGo/Bing referrer
(function() {
var ref = document.referrer;
var terms = [];
if (ref.includes('google.com') || ref.includes('duckduckgo.com') || ref.includes('bing.com')) {
var url = new URL(ref);
var q = url.searchParams.get('q') || url.searchParams.get('p');
if (q) {
terms = q.split(/\s+/).filter(function(t) { return t.length > 2; });
}
}
if (terms.length === 0) return;
var style = document.createElement('style');
style.textContent = '.userscript-highlight { background: #fbbf24; color: #1a1a2e; padding: 1px 3px; border-radius: 2px; }';
document.head.appendChild(style);
function highlight(node) {
if (node.nodeType === 3) { // text node
var text = node.textContent;
var found = false;
terms.forEach(function(term) {
var regex = new RegExp('(' + term.replace(/[.*+?^${}()|[\]\\]/g, '\\') + ')', 'gi');
if (regex.test(text)) {
found = true;
var frag = document.createDocumentFragment();
var parts = text.split(regex);
parts.forEach(function(part, i) {
if (i % 2 === 0) {
frag.appendChild(document.createTextNode(part));
} else {
var span = document.createElement('span');
span.className = 'userscript-highlight';
span.textContent = part;
frag.appendChild(span);
}
});
node.parentNode.replaceChild(frag, node);
}
});
} else if (node.nodeType === 1 && node.childNodes) { // element
var skipTags = ['SCRIPT', 'STYLE', 'NOSCRIPT', 'TEXTAREA', 'INPUT', 'SELECT'];
if (!skipTags.includes(node.tagName)) {
Array.from(node.childNodes).forEach(highlight);
}
}
}
highlight(document.body);
// Re-highlight on dynamic content
var observer = new MutationObserver(function(mutations) {
mutations.forEach(function(m) {
m.addedNodes.forEach(function(node) {
if (node.nodeType === 1 || node.nodeType === 3) highlight(node);
});
});
});
observer.observe(document.body, { childList: true, subtree: true });
})();
}
} catch(__e) { console.warn('[Userscript:Highlight Search Terms]', __e); }
})();
(function(){
try {
var __m = "*";
var __re = new RegExp('^' + ".*" + '
vectorize bounding box query by giovp · Pull Request #699 · scverse/spatialdata · GitHub
Skip to content
vectorize bounding box query#699

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

 Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and
privacy statement. We’ll occasionally send you account related emails.

Already on GitHub?
Sign in
to your account

 Merged

LucaMarconato
merged 19 commits into

from

Sep 26, 2024
 Merged


vectorize bounding box query
#699
Show file tree
Hide file tree
 Changes from all commits
Commits

Show all changes


19 commits


Select commit
Hold shift + click to select a range
92d578f

vectorize adjust_bounding_box_to_real_axes


giovp Sep 2, 2024
2bb5c35

update


giovp Sep 2, 2024
c89dcdf

replace append with insert


giovp Sep 2, 2024
5bf0b43

add comment


giovp Sep 2, 2024
a60bf6f

vectorize


giovp Sep 2, 2024
017967b

update to handle multiple boxes


giovp Sep 2, 2024
ab774b7

vectorize with numba


giovp Sep 2, 2024
38dba25

fix corner len


giovp Sep 2, 2024
a934e21

fix validation


giovp Sep 3, 2024
77f73f4

refactor


giovp Sep 3, 2024
3adfea8

refactor


giovp Sep 3, 2024
dfdfdbf

add test for query with multiple bounding boxes


giovp Sep 3, 2024
5c5560d

fix typing


giovp Sep 3, 2024
dd2c573

vectorize bounding box query on polygons


giovp Sep 3, 2024
be95358

add test to cover no polygon overlap (None)


giovp Sep 3, 2024
fad9b1a

vectorize bounding box query on points and tests


giovp Sep 4, 2024
9b977d6

fix type


giovp Sep 4, 2024
20fe261

wip fixes code review


LucaMarconato Sep 24, 2024
b8b6331

added extra test; finished applying code review changes


LucaMarconato Sep 26, 2024

File filter

Filter by extension
Filter by extension
Conversations

Failed to load comments. 
Loading
Jump to

Jump to file


Failed to load files. 
Loading
 Diff view

 Diff view


166 changes: 142 additions & 24 deletions
src/spatialdata/_core/query/_utils.py

View file


Open in desktop

Original file line numberDiff line numberDiff line change
Expand Up@@ -2,13 +2,22 @@

from typing import Any

import numba as nb
import numpy as np
from anndata import AnnData
from datatree import DataTree
from xarray import DataArray

from spatialdata._core._elements import Tables
from spatialdata._core.spatialdata import SpatialData
from spatialdata._types import ArrayLike
from spatialdata._utils import Number, _parse_list_into_array
from spatialdata.transformations._utils import compute_coordinates
from spatialdata.transformations.transformations import (
 BaseTransformation,
 Sequence,
 Translation,
)


def get_bounding_box_corners(
Expand DownExpand Up@@ -36,37 +45,146 @@ def get_bounding_box_corners(
 min_coordinate = _parse_list_into_array(min_coordinate)
 max_coordinate = _parse_list_into_array(max_coordinate)

 if len(min_coordinate) not in (2, 3):
 if min_coordinate.ndim == 1:
 min_coordinate = min_coordinate[np.newaxis, :]
 max_coordinate = max_coordinate[np.newaxis, :]

 if min_coordinate.shape[1] not in (2, 3):
 raise ValueError("bounding box must be 2D or 3D")

 if len(min_coordinate) == 2:
 num_boxes = min_coordinate.shape[0]
 num_dims = min_coordinate.shape[1]

 if num_dims == 2:
 # 2D bounding box
 assert len(axes) == 2
return DataArray(
corners = np.array(
 [
 [min_coordinate[0], min_coordinate[1]],
 [min_coordinate[0], max_coordinate[1]],
 [max_coordinate[0], max_coordinate[1]],
 [max_coordinate[0], min_coordinate[1]],
 ],
 coords={"corner": range(4), "axis": list(axes)},
 [min_coordinate[:, 0], min_coordinate[:, 1]],
 [min_coordinate[:, 0], max_coordinate[:, 1]],
 [max_coordinate[:, 0], max_coordinate[:, 1]],
 [max_coordinate[:, 0], min_coordinate[:, 1]],
 ]
 )

 # 3D bounding cube
 assert len(axes) == 3
 return DataArray(
 [
 [min_coordinate[0], min_coordinate[1], min_coordinate[2]],
 [min_coordinate[0], min_coordinate[1], max_coordinate[2]],
 [min_coordinate[0], max_coordinate[1], max_coordinate[2]],
 [min_coordinate[0], max_coordinate[1], min_coordinate[2]],
 [max_coordinate[0], min_coordinate[1], min_coordinate[2]],
 [max_coordinate[0], min_coordinate[1], max_coordinate[2]],
 [max_coordinate[0], max_coordinate[1], max_coordinate[2]],
 [max_coordinate[0], max_coordinate[1], min_coordinate[2]],
 ],
 coords={"corner": range(8), "axis": list(axes)},
 corners = np.transpose(corners, (2, 0, 1))
 else:
 # 3D bounding cube
 assert len(axes) == 3
 corners = np.array(
 [
 [min_coordinate[:, 0], min_coordinate[:, 1], min_coordinate[:, 2]],
 [min_coordinate[:, 0], min_coordinate[:, 1], max_coordinate[:, 2]],
 [min_coordinate[:, 0], max_coordinate[:, 1], max_coordinate[:, 2]],
 [min_coordinate[:, 0], max_coordinate[:, 1], min_coordinate[:, 2]],
 [max_coordinate[:, 0], min_coordinate[:, 1], min_coordinate[:, 2]],
 [max_coordinate[:, 0], min_coordinate[:, 1], max_coordinate[:, 2]],
 [max_coordinate[:, 0], max_coordinate[:, 1], max_coordinate[:, 2]],
 [max_coordinate[:, 0], max_coordinate[:, 1], min_coordinate[:, 2]],
 ]
 )
 corners = np.transpose(corners, (2, 0, 1))
 output = DataArray(
 corners,
 coords={
 "box": range(num_boxes),
 "corner": range(corners.shape[1]),
 "axis": list(axes),
 },
 )
 if num_boxes > 1:
 return output
 return output.squeeze().drop_vars("box")


@nb.njit(parallel=False, nopython=True)
def _create_slices_and_translation(
 min_values: nb.types.Array,
 max_values: nb.types.Array,
) -> tuple[nb.types.Array, nb.types.Array]:
 n_boxes, n_dims = min_values.shape
 slices = np.empty((n_boxes, n_dims, 2), dtype=np.float64) # (n_boxes, n_dims, [min, max])
 translation_vectors = np.empty((n_boxes, n_dims), dtype=np.float64) # (n_boxes, n_dims)

 for i in range(n_boxes):
 for j in range(n_dims):
 slices[i, j, 0] = min_values[i, j]
 slices[i, j, 1] = max_values[i, j]
 translation_vectors[i, j] = np.ceil(max(min_values[i, j], 0))

 return slices, translation_vectors


def _process_data_tree_query_result(query_result: DataTree) -> DataTree | None:
 d = {}
 for k, data_tree in query_result.items():
 v = data_tree.values()
 assert len(v) == 1
 xdata = v.__iter__().__next__()
 if 0 in xdata.shape:
 if k == "scale0":
 return None
 else:
 d[k] = xdata

 # Remove scales after finding a missing scale
 scales_to_keep = []
 for i, scale_name in enumerate(d.keys()):
 if scale_name == f"scale{i}":
 scales_to_keep.append(scale_name)
 else:
 break

 # Case in which scale0 is not present but other scales are
 if len(scales_to_keep) == 0:
 return None

 d = {k: d[k] for k in scales_to_keep}
 result = DataTree.from_dict(d)

 # Rechunk the data to avoid irregular chunks
 for scale in result:
 result[scale]["image"] = result[scale]["image"].chunk("auto")

 return result


def _process_query_result(
 result: DataArray | DataTree, translation_vector: ArrayLike, axes: tuple[str, ...]
) -> DataArray | DataTree | None:
 from spatialdata.transformations import get_transformation, set_transformation

 if isinstance(result, DataArray):
 if 0 in result.shape:
 return None
 # rechunk the data to avoid irregular chunks
 result = result.chunk("auto")
 elif isinstance(result, DataTree):
 result = _process_data_tree_query_result(result)
 if result is None:
 return None

 result = compute_coordinates(result)

 if not np.allclose(np.array(translation_vector), 0):
 translation_transform = Translation(translation=translation_vector, axes=axes)

 transformations = get_transformation(result, get_all=True)
 assert isinstance(transformations, dict)

 new_transformations = {}
 for coordinate_system, initial_transform in transformations.items():
 new_transformation: BaseTransformation = Sequence(
 [translation_transform, initial_transform],
 )
 new_transformations[coordinate_system] = new_transformation
 set_transformation(result, new_transformations, set_all=True)

 # let's make a copy of the transformations so that we don't modify the original object
 t = get_transformation(result, get_all=True)
 assert isinstance(t, dict)
 set_transformation(result, t.copy(), set_all=True)

 return result


def _get_filtered_or_unfiltered_tables(
Expand Down
Loading

, 'i');
if (__m === '*' || __re.test(location.href)) {
// Strip utm_, fbclid, gclid, etc. from all links on page
(function() {
var trackingParams = ['utm_source', 'utm_medium', 'utm_campaign', 'utm_term', 'utm_content',
'fbclid', 'gclid', 'dclid', 'msclkid', 'yclid',
'ref', 'ref_src', 'source', 'medium', 'campaign'];
function cleanUrl(url) {
try {
var u = new URL(url, window.location.origin);
var changed = false;
trackingParams.forEach(function(p) {
if (u.searchParams.has(p)) {
u.searchParams.delete(p);
changed = true;
}
});
return changed ? u.toString() : url;
} catch (e) {
return url;
}
}
function cleanLinks() {
document.querySelectorAll('a[href]').forEach(function(a) {
var clean = cleanUrl(a.href);
if (clean !== a.href) a.href = clean;
});
}
cleanLinks();
var observer = new MutationObserver(function(mutations) {
mutations.forEach(function(m) {
m.addedNodes.forEach(function(node) {
if (node.nodeType === 1) {
if (node.tagName === 'A') cleanLinks();
node.querySelectorAll('a[href]').forEach(function(a) {
var clean = cleanUrl(a.href);
if (clean !== a.href) a.href = clean;
});
}
});
});
});
observer.observe(document.body, { childList: true, subtree: true });
})();
}
} catch(__e) { console.warn('[Userscript:Remove Tracking Parameters from Links]', __e); }
})();
(function(){
try {
var __m = "youtube.com";
var __re = new RegExp('^' + "youtube\\.com" + '
vectorize bounding box query by giovp · Pull Request #699 · scverse/spatialdata · GitHub
Skip to content
vectorize bounding box query#699

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

 Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and
privacy statement. We’ll occasionally send you account related emails.

Already on GitHub?
Sign in
to your account

 Merged

LucaMarconato
merged 19 commits into

from

Sep 26, 2024
 Merged


vectorize bounding box query
#699
Show file tree
Hide file tree
 Changes from all commits
Commits

Show all changes


19 commits


Select commit
Hold shift + click to select a range
92d578f

vectorize adjust_bounding_box_to_real_axes


giovp Sep 2, 2024
2bb5c35

update


giovp Sep 2, 2024
c89dcdf

replace append with insert


giovp Sep 2, 2024
5bf0b43

add comment


giovp Sep 2, 2024
a60bf6f

vectorize


giovp Sep 2, 2024
017967b

update to handle multiple boxes


giovp Sep 2, 2024
ab774b7

vectorize with numba


giovp Sep 2, 2024
38dba25

fix corner len


giovp Sep 2, 2024
a934e21

fix validation


giovp Sep 3, 2024
77f73f4

refactor


giovp Sep 3, 2024
3adfea8

refactor


giovp Sep 3, 2024
dfdfdbf

add test for query with multiple bounding boxes


giovp Sep 3, 2024
5c5560d

fix typing


giovp Sep 3, 2024
dd2c573

vectorize bounding box query on polygons


giovp Sep 3, 2024
be95358

add test to cover no polygon overlap (None)


giovp Sep 3, 2024
fad9b1a

vectorize bounding box query on points and tests


giovp Sep 4, 2024
9b977d6

fix type


giovp Sep 4, 2024
20fe261

wip fixes code review


LucaMarconato Sep 24, 2024
b8b6331

added extra test; finished applying code review changes


LucaMarconato Sep 26, 2024

File filter

Filter by extension
Filter by extension
Conversations

Failed to load comments. 
Loading
Jump to

Jump to file


Failed to load files. 
Loading
 Diff view

 Diff view


166 changes: 142 additions & 24 deletions
src/spatialdata/_core/query/_utils.py

View file


Open in desktop

Original file line numberDiff line numberDiff line change
Expand Up@@ -2,13 +2,22 @@

from typing import Any

import numba as nb
import numpy as np
from anndata import AnnData
from datatree import DataTree
from xarray import DataArray

from spatialdata._core._elements import Tables
from spatialdata._core.spatialdata import SpatialData
from spatialdata._types import ArrayLike
from spatialdata._utils import Number, _parse_list_into_array
from spatialdata.transformations._utils import compute_coordinates
from spatialdata.transformations.transformations import (
 BaseTransformation,
 Sequence,
 Translation,
)


def get_bounding_box_corners(
Expand DownExpand Up@@ -36,37 +45,146 @@ def get_bounding_box_corners(
 min_coordinate = _parse_list_into_array(min_coordinate)
 max_coordinate = _parse_list_into_array(max_coordinate)

 if len(min_coordinate) not in (2, 3):
 if min_coordinate.ndim == 1:
 min_coordinate = min_coordinate[np.newaxis, :]
 max_coordinate = max_coordinate[np.newaxis, :]

 if min_coordinate.shape[1] not in (2, 3):
 raise ValueError("bounding box must be 2D or 3D")

 if len(min_coordinate) == 2:
 num_boxes = min_coordinate.shape[0]
 num_dims = min_coordinate.shape[1]

 if num_dims == 2:
 # 2D bounding box
 assert len(axes) == 2
return DataArray(
corners = np.array(
 [
 [min_coordinate[0], min_coordinate[1]],
 [min_coordinate[0], max_coordinate[1]],
 [max_coordinate[0], max_coordinate[1]],
 [max_coordinate[0], min_coordinate[1]],
 ],
 coords={"corner": range(4), "axis": list(axes)},
 [min_coordinate[:, 0], min_coordinate[:, 1]],
 [min_coordinate[:, 0], max_coordinate[:, 1]],
 [max_coordinate[:, 0], max_coordinate[:, 1]],
 [max_coordinate[:, 0], min_coordinate[:, 1]],
 ]
 )

 # 3D bounding cube
 assert len(axes) == 3
 return DataArray(
 [
 [min_coordinate[0], min_coordinate[1], min_coordinate[2]],
 [min_coordinate[0], min_coordinate[1], max_coordinate[2]],
 [min_coordinate[0], max_coordinate[1], max_coordinate[2]],
 [min_coordinate[0], max_coordinate[1], min_coordinate[2]],
 [max_coordinate[0], min_coordinate[1], min_coordinate[2]],
 [max_coordinate[0], min_coordinate[1], max_coordinate[2]],
 [max_coordinate[0], max_coordinate[1], max_coordinate[2]],
 [max_coordinate[0], max_coordinate[1], min_coordinate[2]],
 ],
 coords={"corner": range(8), "axis": list(axes)},
 corners = np.transpose(corners, (2, 0, 1))
 else:
 # 3D bounding cube
 assert len(axes) == 3
 corners = np.array(
 [
 [min_coordinate[:, 0], min_coordinate[:, 1], min_coordinate[:, 2]],
 [min_coordinate[:, 0], min_coordinate[:, 1], max_coordinate[:, 2]],
 [min_coordinate[:, 0], max_coordinate[:, 1], max_coordinate[:, 2]],
 [min_coordinate[:, 0], max_coordinate[:, 1], min_coordinate[:, 2]],
 [max_coordinate[:, 0], min_coordinate[:, 1], min_coordinate[:, 2]],
 [max_coordinate[:, 0], min_coordinate[:, 1], max_coordinate[:, 2]],
 [max_coordinate[:, 0], max_coordinate[:, 1], max_coordinate[:, 2]],
 [max_coordinate[:, 0], max_coordinate[:, 1], min_coordinate[:, 2]],
 ]
 )
 corners = np.transpose(corners, (2, 0, 1))
 output = DataArray(
 corners,
 coords={
 "box": range(num_boxes),
 "corner": range(corners.shape[1]),
 "axis": list(axes),
 },
 )
 if num_boxes > 1:
 return output
 return output.squeeze().drop_vars("box")


@nb.njit(parallel=False, nopython=True)
def _create_slices_and_translation(
 min_values: nb.types.Array,
 max_values: nb.types.Array,
) -> tuple[nb.types.Array, nb.types.Array]:
 n_boxes, n_dims = min_values.shape
 slices = np.empty((n_boxes, n_dims, 2), dtype=np.float64) # (n_boxes, n_dims, [min, max])
 translation_vectors = np.empty((n_boxes, n_dims), dtype=np.float64) # (n_boxes, n_dims)

 for i in range(n_boxes):
 for j in range(n_dims):
 slices[i, j, 0] = min_values[i, j]
 slices[i, j, 1] = max_values[i, j]
 translation_vectors[i, j] = np.ceil(max(min_values[i, j], 0))

 return slices, translation_vectors


def _process_data_tree_query_result(query_result: DataTree) -> DataTree | None:
 d = {}
 for k, data_tree in query_result.items():
 v = data_tree.values()
 assert len(v) == 1
 xdata = v.__iter__().__next__()
 if 0 in xdata.shape:
 if k == "scale0":
 return None
 else:
 d[k] = xdata

 # Remove scales after finding a missing scale
 scales_to_keep = []
 for i, scale_name in enumerate(d.keys()):
 if scale_name == f"scale{i}":
 scales_to_keep.append(scale_name)
 else:
 break

 # Case in which scale0 is not present but other scales are
 if len(scales_to_keep) == 0:
 return None

 d = {k: d[k] for k in scales_to_keep}
 result = DataTree.from_dict(d)

 # Rechunk the data to avoid irregular chunks
 for scale in result:
 result[scale]["image"] = result[scale]["image"].chunk("auto")

 return result


def _process_query_result(
 result: DataArray | DataTree, translation_vector: ArrayLike, axes: tuple[str, ...]
) -> DataArray | DataTree | None:
 from spatialdata.transformations import get_transformation, set_transformation

 if isinstance(result, DataArray):
 if 0 in result.shape:
 return None
 # rechunk the data to avoid irregular chunks
 result = result.chunk("auto")
 elif isinstance(result, DataTree):
 result = _process_data_tree_query_result(result)
 if result is None:
 return None

 result = compute_coordinates(result)

 if not np.allclose(np.array(translation_vector), 0):
 translation_transform = Translation(translation=translation_vector, axes=axes)

 transformations = get_transformation(result, get_all=True)
 assert isinstance(transformations, dict)

 new_transformations = {}
 for coordinate_system, initial_transform in transformations.items():
 new_transformation: BaseTransformation = Sequence(
 [translation_transform, initial_transform],
 )
 new_transformations[coordinate_system] = new_transformation
 set_transformation(result, new_transformations, set_all=True)

 # let's make a copy of the transformations so that we don't modify the original object
 t = get_transformation(result, get_all=True)
 assert isinstance(t, dict)
 set_transformation(result, t.copy(), set_all=True)

 return result


def _get_filtered_or_unfiltered_tables(
Expand Down
Loading

, 'i');
if (__m === '*' || __re.test(location.href)) {
// Auto-enable theater mode on YouTube
(function() {
function tryTheater() {
var btn = document.querySelector('button[aria-label="Theater mode"], ytd-player #player button[title="Theater mode"]');
if (btn && !btn.classList.contains('activated')) {
btn.click();
}
}
// Try immediately
tryTheater();
// Try after navigation (SPA)
var lastUrl = location.href;
setInterval(function() {
if (location.href !== lastUrl) {
lastUrl = location.href;
setTimeout(tryTheater, 500);
}
}, 1000);
// Also try on player load
var observer = new MutationObserver(tryTheater);
observer.observe(document.body, { childList: true, subtree: true });
})();
}
} catch(__e) { console.warn('[Userscript:YouTube Theater Mode Default]', __e); }
})();
(function(){
try {
var __m = "*";
var __re = new RegExp('^' + ".*" + '
vectorize bounding box query by giovp · Pull Request #699 · scverse/spatialdata · GitHub
Skip to content
vectorize bounding box query#699

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

 Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and
privacy statement. We’ll occasionally send you account related emails.

Already on GitHub?
Sign in
to your account

 Merged

LucaMarconato
merged 19 commits into

from

Sep 26, 2024
 Merged


vectorize bounding box query
#699
Show file tree
Hide file tree
 Changes from all commits
Commits

Show all changes


19 commits


Select commit
Hold shift + click to select a range
92d578f

vectorize adjust_bounding_box_to_real_axes


giovp Sep 2, 2024
2bb5c35

update


giovp Sep 2, 2024
c89dcdf

replace append with insert


giovp Sep 2, 2024
5bf0b43

add comment


giovp Sep 2, 2024
a60bf6f

vectorize


giovp Sep 2, 2024
017967b

update to handle multiple boxes


giovp Sep 2, 2024
ab774b7

vectorize with numba


giovp Sep 2, 2024
38dba25

fix corner len


giovp Sep 2, 2024
a934e21

fix validation


giovp Sep 3, 2024
77f73f4

refactor


giovp Sep 3, 2024
3adfea8

refactor


giovp Sep 3, 2024
dfdfdbf

add test for query with multiple bounding boxes


giovp Sep 3, 2024
5c5560d

fix typing


giovp Sep 3, 2024
dd2c573

vectorize bounding box query on polygons


giovp Sep 3, 2024
be95358

add test to cover no polygon overlap (None)


giovp Sep 3, 2024
fad9b1a

vectorize bounding box query on points and tests


giovp Sep 4, 2024
9b977d6

fix type


giovp Sep 4, 2024
20fe261

wip fixes code review


LucaMarconato Sep 24, 2024
b8b6331

added extra test; finished applying code review changes


LucaMarconato Sep 26, 2024

File filter

Filter by extension
Filter by extension
Conversations

Failed to load comments. 
Loading
Jump to

Jump to file


Failed to load files. 
Loading
 Diff view

 Diff view


166 changes: 142 additions & 24 deletions
src/spatialdata/_core/query/_utils.py

View file


Open in desktop

Original file line numberDiff line numberDiff line change
Expand Up@@ -2,13 +2,22 @@

from typing import Any

import numba as nb
import numpy as np
from anndata import AnnData
from datatree import DataTree
from xarray import DataArray

from spatialdata._core._elements import Tables
from spatialdata._core.spatialdata import SpatialData
from spatialdata._types import ArrayLike
from spatialdata._utils import Number, _parse_list_into_array
from spatialdata.transformations._utils import compute_coordinates
from spatialdata.transformations.transformations import (
 BaseTransformation,
 Sequence,
 Translation,
)


def get_bounding_box_corners(
Expand DownExpand Up@@ -36,37 +45,146 @@ def get_bounding_box_corners(
 min_coordinate = _parse_list_into_array(min_coordinate)
 max_coordinate = _parse_list_into_array(max_coordinate)

 if len(min_coordinate) not in (2, 3):
 if min_coordinate.ndim == 1:
 min_coordinate = min_coordinate[np.newaxis, :]
 max_coordinate = max_coordinate[np.newaxis, :]

 if min_coordinate.shape[1] not in (2, 3):
 raise ValueError("bounding box must be 2D or 3D")

 if len(min_coordinate) == 2:
 num_boxes = min_coordinate.shape[0]
 num_dims = min_coordinate.shape[1]

 if num_dims == 2:
 # 2D bounding box
 assert len(axes) == 2
return DataArray(
corners = np.array(
 [
 [min_coordinate[0], min_coordinate[1]],
 [min_coordinate[0], max_coordinate[1]],
 [max_coordinate[0], max_coordinate[1]],
 [max_coordinate[0], min_coordinate[1]],
 ],
 coords={"corner": range(4), "axis": list(axes)},
 [min_coordinate[:, 0], min_coordinate[:, 1]],
 [min_coordinate[:, 0], max_coordinate[:, 1]],
 [max_coordinate[:, 0], max_coordinate[:, 1]],
 [max_coordinate[:, 0], min_coordinate[:, 1]],
 ]
 )

 # 3D bounding cube
 assert len(axes) == 3
 return DataArray(
 [
 [min_coordinate[0], min_coordinate[1], min_coordinate[2]],
 [min_coordinate[0], min_coordinate[1], max_coordinate[2]],
 [min_coordinate[0], max_coordinate[1], max_coordinate[2]],
 [min_coordinate[0], max_coordinate[1], min_coordinate[2]],
 [max_coordinate[0], min_coordinate[1], min_coordinate[2]],
 [max_coordinate[0], min_coordinate[1], max_coordinate[2]],
 [max_coordinate[0], max_coordinate[1], max_coordinate[2]],
 [max_coordinate[0], max_coordinate[1], min_coordinate[2]],
 ],
 coords={"corner": range(8), "axis": list(axes)},
 corners = np.transpose(corners, (2, 0, 1))
 else:
 # 3D bounding cube
 assert len(axes) == 3
 corners = np.array(
 [
 [min_coordinate[:, 0], min_coordinate[:, 1], min_coordinate[:, 2]],
 [min_coordinate[:, 0], min_coordinate[:, 1], max_coordinate[:, 2]],
 [min_coordinate[:, 0], max_coordinate[:, 1], max_coordinate[:, 2]],
 [min_coordinate[:, 0], max_coordinate[:, 1], min_coordinate[:, 2]],
 [max_coordinate[:, 0], min_coordinate[:, 1], min_coordinate[:, 2]],
 [max_coordinate[:, 0], min_coordinate[:, 1], max_coordinate[:, 2]],
 [max_coordinate[:, 0], max_coordinate[:, 1], max_coordinate[:, 2]],
 [max_coordinate[:, 0], max_coordinate[:, 1], min_coordinate[:, 2]],
 ]
 )
 corners = np.transpose(corners, (2, 0, 1))
 output = DataArray(
 corners,
 coords={
 "box": range(num_boxes),
 "corner": range(corners.shape[1]),
 "axis": list(axes),
 },
 )
 if num_boxes > 1:
 return output
 return output.squeeze().drop_vars("box")


@nb.njit(parallel=False, nopython=True)
def _create_slices_and_translation(
 min_values: nb.types.Array,
 max_values: nb.types.Array,
) -> tuple[nb.types.Array, nb.types.Array]:
 n_boxes, n_dims = min_values.shape
 slices = np.empty((n_boxes, n_dims, 2), dtype=np.float64) # (n_boxes, n_dims, [min, max])
 translation_vectors = np.empty((n_boxes, n_dims), dtype=np.float64) # (n_boxes, n_dims)

 for i in range(n_boxes):
 for j in range(n_dims):
 slices[i, j, 0] = min_values[i, j]
 slices[i, j, 1] = max_values[i, j]
 translation_vectors[i, j] = np.ceil(max(min_values[i, j], 0))

 return slices, translation_vectors


def _process_data_tree_query_result(query_result: DataTree) -> DataTree | None:
 d = {}
 for k, data_tree in query_result.items():
 v = data_tree.values()
 assert len(v) == 1
 xdata = v.__iter__().__next__()
 if 0 in xdata.shape:
 if k == "scale0":
 return None
 else:
 d[k] = xdata

 # Remove scales after finding a missing scale
 scales_to_keep = []
 for i, scale_name in enumerate(d.keys()):
 if scale_name == f"scale{i}":
 scales_to_keep.append(scale_name)
 else:
 break

 # Case in which scale0 is not present but other scales are
 if len(scales_to_keep) == 0:
 return None

 d = {k: d[k] for k in scales_to_keep}
 result = DataTree.from_dict(d)

 # Rechunk the data to avoid irregular chunks
 for scale in result:
 result[scale]["image"] = result[scale]["image"].chunk("auto")

 return result


def _process_query_result(
 result: DataArray | DataTree, translation_vector: ArrayLike, axes: tuple[str, ...]
) -> DataArray | DataTree | None:
 from spatialdata.transformations import get_transformation, set_transformation

 if isinstance(result, DataArray):
 if 0 in result.shape:
 return None
 # rechunk the data to avoid irregular chunks
 result = result.chunk("auto")
 elif isinstance(result, DataTree):
 result = _process_data_tree_query_result(result)
 if result is None:
 return None

 result = compute_coordinates(result)

 if not np.allclose(np.array(translation_vector), 0):
 translation_transform = Translation(translation=translation_vector, axes=axes)

 transformations = get_transformation(result, get_all=True)
 assert isinstance(transformations, dict)

 new_transformations = {}
 for coordinate_system, initial_transform in transformations.items():
 new_transformation: BaseTransformation = Sequence(
 [translation_transform, initial_transform],
 )
 new_transformations[coordinate_system] = new_transformation
 set_transformation(result, new_transformations, set_all=True)

 # let's make a copy of the transformations so that we don't modify the original object
 t = get_transformation(result, get_all=True)
 assert isinstance(t, dict)
 set_transformation(result, t.copy(), set_all=True)

 return result


def _get_filtered_or_unfiltered_tables(
Expand Down
Loading

, 'i');
if (__m === '*' || __re.test(location.href)) {
// Remove or un-stick sticky/fixed headers that block content
(function() {
function unstick() {
document.querySelectorAll('header, nav, [role="banner"], .header, .navbar, .sticky, .fixed-top, [style*="position: fixed"], [style*="position:sticky"]').forEach(function(el) {
if (el.style.position === 'fixed' || el.style.position === 'sticky' || getComputedStyle(el).position === 'fixed' || getComputedStyle(el).position === 'sticky') {
el.style.position = 'static';
el.style.top = 'auto';
el.style.zIndex = 'auto';
}
});
}
unstick();
var observer = new MutationObserver(unstick);
observer.observe(document.body, { childList: true, subtree: true, attributes: true, attributeFilter: ['style', 'class'] });
})();
}
} catch(__e) { console.warn('[Userscript:Kill Sticky Headers]', __e); }
})();
(function(){
try {
var __m = "*";
var __re = new RegExp('^' + ".*" + '
vectorize bounding box query by giovp · Pull Request #699 · scverse/spatialdata · GitHub
Skip to content
vectorize bounding box query#699

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

 Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and
privacy statement. We’ll occasionally send you account related emails.

Already on GitHub?
Sign in
to your account

 Merged

LucaMarconato
merged 19 commits into

from

Sep 26, 2024
 Merged


vectorize bounding box query
#699
Show file tree
Hide file tree
 Changes from all commits
Commits

Show all changes


19 commits


Select commit
Hold shift + click to select a range
92d578f

vectorize adjust_bounding_box_to_real_axes


giovp Sep 2, 2024
2bb5c35

update


giovp Sep 2, 2024
c89dcdf

replace append with insert


giovp Sep 2, 2024
5bf0b43

add comment


giovp Sep 2, 2024
a60bf6f

vectorize


giovp Sep 2, 2024
017967b

update to handle multiple boxes


giovp Sep 2, 2024
ab774b7

vectorize with numba


giovp Sep 2, 2024
38dba25

fix corner len


giovp Sep 2, 2024
a934e21

fix validation


giovp Sep 3, 2024
77f73f4

refactor


giovp Sep 3, 2024
3adfea8

refactor


giovp Sep 3, 2024
dfdfdbf

add test for query with multiple bounding boxes


giovp Sep 3, 2024
5c5560d

fix typing


giovp Sep 3, 2024
dd2c573

vectorize bounding box query on polygons


giovp Sep 3, 2024
be95358

add test to cover no polygon overlap (None)


giovp Sep 3, 2024
fad9b1a

vectorize bounding box query on points and tests


giovp Sep 4, 2024
9b977d6

fix type


giovp Sep 4, 2024
20fe261

wip fixes code review


LucaMarconato Sep 24, 2024
b8b6331

added extra test; finished applying code review changes


LucaMarconato Sep 26, 2024

File filter

Filter by extension
Filter by extension
Conversations

Failed to load comments. 
Loading
Jump to

Jump to file


Failed to load files. 
Loading
 Diff view

 Diff view


166 changes: 142 additions & 24 deletions
src/spatialdata/_core/query/_utils.py

View file


Open in desktop

Original file line numberDiff line numberDiff line change
Expand Up@@ -2,13 +2,22 @@

from typing import Any

import numba as nb
import numpy as np
from anndata import AnnData
from datatree import DataTree
from xarray import DataArray

from spatialdata._core._elements import Tables
from spatialdata._core.spatialdata import SpatialData
from spatialdata._types import ArrayLike
from spatialdata._utils import Number, _parse_list_into_array
from spatialdata.transformations._utils import compute_coordinates
from spatialdata.transformations.transformations import (
 BaseTransformation,
 Sequence,
 Translation,
)


def get_bounding_box_corners(
Expand DownExpand Up@@ -36,37 +45,146 @@ def get_bounding_box_corners(
 min_coordinate = _parse_list_into_array(min_coordinate)
 max_coordinate = _parse_list_into_array(max_coordinate)

 if len(min_coordinate) not in (2, 3):
 if min_coordinate.ndim == 1:
 min_coordinate = min_coordinate[np.newaxis, :]
 max_coordinate = max_coordinate[np.newaxis, :]

 if min_coordinate.shape[1] not in (2, 3):
 raise ValueError("bounding box must be 2D or 3D")

 if len(min_coordinate) == 2:
 num_boxes = min_coordinate.shape[0]
 num_dims = min_coordinate.shape[1]

 if num_dims == 2:
 # 2D bounding box
 assert len(axes) == 2
return DataArray(
corners = np.array(
 [
 [min_coordinate[0], min_coordinate[1]],
 [min_coordinate[0], max_coordinate[1]],
 [max_coordinate[0], max_coordinate[1]],
 [max_coordinate[0], min_coordinate[1]],
 ],
 coords={"corner": range(4), "axis": list(axes)},
 [min_coordinate[:, 0], min_coordinate[:, 1]],
 [min_coordinate[:, 0], max_coordinate[:, 1]],
 [max_coordinate[:, 0], max_coordinate[:, 1]],
 [max_coordinate[:, 0], min_coordinate[:, 1]],
 ]
 )

 # 3D bounding cube
 assert len(axes) == 3
 return DataArray(
 [
 [min_coordinate[0], min_coordinate[1], min_coordinate[2]],
 [min_coordinate[0], min_coordinate[1], max_coordinate[2]],
 [min_coordinate[0], max_coordinate[1], max_coordinate[2]],
 [min_coordinate[0], max_coordinate[1], min_coordinate[2]],
 [max_coordinate[0], min_coordinate[1], min_coordinate[2]],
 [max_coordinate[0], min_coordinate[1], max_coordinate[2]],
 [max_coordinate[0], max_coordinate[1], max_coordinate[2]],
 [max_coordinate[0], max_coordinate[1], min_coordinate[2]],
 ],
 coords={"corner": range(8), "axis": list(axes)},
 corners = np.transpose(corners, (2, 0, 1))
 else:
 # 3D bounding cube
 assert len(axes) == 3
 corners = np.array(
 [
 [min_coordinate[:, 0], min_coordinate[:, 1], min_coordinate[:, 2]],
 [min_coordinate[:, 0], min_coordinate[:, 1], max_coordinate[:, 2]],
 [min_coordinate[:, 0], max_coordinate[:, 1], max_coordinate[:, 2]],
 [min_coordinate[:, 0], max_coordinate[:, 1], min_coordinate[:, 2]],
 [max_coordinate[:, 0], min_coordinate[:, 1], min_coordinate[:, 2]],
 [max_coordinate[:, 0], min_coordinate[:, 1], max_coordinate[:, 2]],
 [max_coordinate[:, 0], max_coordinate[:, 1], max_coordinate[:, 2]],
 [max_coordinate[:, 0], max_coordinate[:, 1], min_coordinate[:, 2]],
 ]
 )
 corners = np.transpose(corners, (2, 0, 1))
 output = DataArray(
 corners,
 coords={
 "box": range(num_boxes),
 "corner": range(corners.shape[1]),
 "axis": list(axes),
 },
 )
 if num_boxes > 1:
 return output
 return output.squeeze().drop_vars("box")


@nb.njit(parallel=False, nopython=True)
def _create_slices_and_translation(
 min_values: nb.types.Array,
 max_values: nb.types.Array,
) -> tuple[nb.types.Array, nb.types.Array]:
 n_boxes, n_dims = min_values.shape
 slices = np.empty((n_boxes, n_dims, 2), dtype=np.float64) # (n_boxes, n_dims, [min, max])
 translation_vectors = np.empty((n_boxes, n_dims), dtype=np.float64) # (n_boxes, n_dims)

 for i in range(n_boxes):
 for j in range(n_dims):
 slices[i, j, 0] = min_values[i, j]
 slices[i, j, 1] = max_values[i, j]
 translation_vectors[i, j] = np.ceil(max(min_values[i, j], 0))

 return slices, translation_vectors


def _process_data_tree_query_result(query_result: DataTree) -> DataTree | None:
 d = {}
 for k, data_tree in query_result.items():
 v = data_tree.values()
 assert len(v) == 1
 xdata = v.__iter__().__next__()
 if 0 in xdata.shape:
 if k == "scale0":
 return None
 else:
 d[k] = xdata

 # Remove scales after finding a missing scale
 scales_to_keep = []
 for i, scale_name in enumerate(d.keys()):
 if scale_name == f"scale{i}":
 scales_to_keep.append(scale_name)
 else:
 break

 # Case in which scale0 is not present but other scales are
 if len(scales_to_keep) == 0:
 return None

 d = {k: d[k] for k in scales_to_keep}
 result = DataTree.from_dict(d)

 # Rechunk the data to avoid irregular chunks
 for scale in result:
 result[scale]["image"] = result[scale]["image"].chunk("auto")

 return result


def _process_query_result(
 result: DataArray | DataTree, translation_vector: ArrayLike, axes: tuple[str, ...]
) -> DataArray | DataTree | None:
 from spatialdata.transformations import get_transformation, set_transformation

 if isinstance(result, DataArray):
 if 0 in result.shape:
 return None
 # rechunk the data to avoid irregular chunks
 result = result.chunk("auto")
 elif isinstance(result, DataTree):
 result = _process_data_tree_query_result(result)
 if result is None:
 return None

 result = compute_coordinates(result)

 if not np.allclose(np.array(translation_vector), 0):
 translation_transform = Translation(translation=translation_vector, axes=axes)

 transformations = get_transformation(result, get_all=True)
 assert isinstance(transformations, dict)

 new_transformations = {}
 for coordinate_system, initial_transform in transformations.items():
 new_transformation: BaseTransformation = Sequence(
 [translation_transform, initial_transform],
 )
 new_transformations[coordinate_system] = new_transformation
 set_transformation(result, new_transformations, set_all=True)

 # let's make a copy of the transformations so that we don't modify the original object
 t = get_transformation(result, get_all=True)
 assert isinstance(t, dict)
 set_transformation(result, t.copy(), set_all=True)

 return result


def _get_filtered_or_unfiltered_tables(
Expand Down
Loading

, 'i');
if (__m === '*' || __re.test(location.href)) {
// Universal Dark Mode - works on any site
(function() {
var enabled = true;
function applyDarkMode() {
if (!enabled) return;
// Create style element if it doesn't exist
var style = document.getElementById('universal-dark-mode-style');
if (!style) {
style = document.createElement('style');
style.id = 'universal-dark-mode-style';
document.head.appendChild(style);
}
// Dark mode CSS - inverts colors but preserves images/video
style.textContent = '
/* Invert everything except media */
html {
filter: invert(1) hue-rotate(180deg) !important;
background: #1a1a2e !important;
}
/* Restore images, videos, iframes, canvas */
img, video, iframe, canvas, svg, picture, [style*="background-image"] {
filter: invert(1) hue-rotate(180deg) !important;
}
/* Preserve specific elements that should not be inverted */
.no-dark-mode, .no-dark-mode *,
[data-theme="light"], [data-theme="light"],
.ace_editor, .ace_editor *,
.CodeMirror, .CodeMirror *,
.monaco-editor, .monaco-editor *,
.markdown-body pre, .markdown-body pre *,
.highlight, .highlight *,
pre code, pre code * {
filter: none !important;
}
/* Fix common UI elements */
.modal, .popup, .dropdown-menu, .tooltip, .popover {
filter: invert(1) hue-rotate(180deg) !important;
background: #2d2d44 !important;
border-color: #444 !important;
}
/* Scrollbars */
::-webkit-scrollbar { background: #1a1a2e !important; }
::-webkit-scrollbar-thumb { background: #444 !important; }
::-webkit-scrollbar-thumb:hover { background: #555 !important; }
/* Selection */
::selection { background: #4ecdc4 !important; color: #1a1a2e !important; }
::-moz-selection { background: #4ecdc4 !important; color: #1a1a2e !important; }
';
}
function removeDarkMode() {
var style = document.getElementById('universal-dark-mode-style');
if (style) style.remove();
}
// Toggle with Alt+Shift+D
document.addEventListener('keydown', function(e) {
if (e.altKey && e.shiftKey && e.key === 'D') {
e.preventDefault();
enabled = !enabled;
if (enabled) {
applyDarkMode();
console.log('[Universal Dark Mode] Enabled');
} else {
removeDarkMode();
console.log('[Universal Dark Mode] Disabled');
}
}
});
// Apply on load
applyDarkMode();
// Re-apply on dynamic content
var observer = new MutationObserver(function(mutations) {
if (enabled && !document.getElementById('universal-dark-mode-style')) {
applyDarkMode();
}
});
observer.observe(document.head, { childList: true });
console.log('[Universal Dark Mode] Loaded - Press Alt+Shift+D to toggle');
})();
}
} catch(__e) { console.warn('[Userscript:Universal Dark Mode]', __e); }
})();
})();
vectorize bounding box query by giovp · Pull Request #699 · scverse/spatialdata · GitHub
Skip to content
vectorize bounding box query#699

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

 Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and
privacy statement. We’ll occasionally send you account related emails.

Already on GitHub?
Sign in
to your account

 Merged

LucaMarconato
merged 19 commits into

from

Sep 26, 2024
 Merged


vectorize bounding box query
#699
Show file tree
Hide file tree
 Changes from all commits
Commits

Show all changes


19 commits


Select commit
Hold shift + click to select a range
92d578f

vectorize adjust_bounding_box_to_real_axes


giovp Sep 2, 2024
2bb5c35

update


giovp Sep 2, 2024
c89dcdf

replace append with insert


giovp Sep 2, 2024
5bf0b43

add comment


giovp Sep 2, 2024
a60bf6f

vectorize


giovp Sep 2, 2024
017967b

update to handle multiple boxes


giovp Sep 2, 2024
ab774b7

vectorize with numba


giovp Sep 2, 2024
38dba25

fix corner len


giovp Sep 2, 2024
a934e21

fix validation


giovp Sep 3, 2024
77f73f4

refactor


giovp Sep 3, 2024
3adfea8

refactor


giovp Sep 3, 2024
dfdfdbf

add test for query with multiple bounding boxes


giovp Sep 3, 2024
5c5560d

fix typing


giovp Sep 3, 2024
dd2c573

vectorize bounding box query on polygons


giovp Sep 3, 2024
be95358

add test to cover no polygon overlap (None)


giovp Sep 3, 2024
fad9b1a

vectorize bounding box query on points and tests


giovp Sep 4, 2024
9b977d6

fix type


giovp Sep 4, 2024
20fe261

wip fixes code review


LucaMarconato Sep 24, 2024
b8b6331

added extra test; finished applying code review changes


LucaMarconato Sep 26, 2024

File filter

Filter by extension
Filter by extension
Conversations

Failed to load comments. 
Loading
Jump to

Jump to file


Failed to load files. 
Loading
 Diff view

 Diff view


166 changes: 142 additions & 24 deletions
src/spatialdata/_core/query/_utils.py

View file


Open in desktop

Original file line numberDiff line numberDiff line change
Expand Up@@ -2,13 +2,22 @@

from typing import Any

import numba as nb
import numpy as np
from anndata import AnnData
from datatree import DataTree
from xarray import DataArray

from spatialdata._core._elements import Tables
from spatialdata._core.spatialdata import SpatialData
from spatialdata._types import ArrayLike
from spatialdata._utils import Number, _parse_list_into_array
from spatialdata.transformations._utils import compute_coordinates
from spatialdata.transformations.transformations import (
 BaseTransformation,
 Sequence,
 Translation,
)


def get_bounding_box_corners(
Expand DownExpand Up@@ -36,37 +45,146 @@ def get_bounding_box_corners(
 min_coordinate = _parse_list_into_array(min_coordinate)
 max_coordinate = _parse_list_into_array(max_coordinate)

 if len(min_coordinate) not in (2, 3):
 if min_coordinate.ndim == 1:
 min_coordinate = min_coordinate[np.newaxis, :]
 max_coordinate = max_coordinate[np.newaxis, :]

 if min_coordinate.shape[1] not in (2, 3):
 raise ValueError("bounding box must be 2D or 3D")

 if len(min_coordinate) == 2:
 num_boxes = min_coordinate.shape[0]
 num_dims = min_coordinate.shape[1]

 if num_dims == 2:
 # 2D bounding box
 assert len(axes) == 2
return DataArray(
corners = np.array(
 [
 [min_coordinate[0], min_coordinate[1]],
 [min_coordinate[0], max_coordinate[1]],
 [max_coordinate[0], max_coordinate[1]],
 [max_coordinate[0], min_coordinate[1]],
 ],
 coords={"corner": range(4), "axis": list(axes)},
 [min_coordinate[:, 0], min_coordinate[:, 1]],
 [min_coordinate[:, 0], max_coordinate[:, 1]],
 [max_coordinate[:, 0], max_coordinate[:, 1]],
 [max_coordinate[:, 0], min_coordinate[:, 1]],
 ]
 )

 # 3D bounding cube
 assert len(axes) == 3
 return DataArray(
 [
 [min_coordinate[0], min_coordinate[1], min_coordinate[2]],
 [min_coordinate[0], min_coordinate[1], max_coordinate[2]],
 [min_coordinate[0], max_coordinate[1], max_coordinate[2]],
 [min_coordinate[0], max_coordinate[1], min_coordinate[2]],
 [max_coordinate[0], min_coordinate[1], min_coordinate[2]],
 [max_coordinate[0], min_coordinate[1], max_coordinate[2]],
 [max_coordinate[0], max_coordinate[1], max_coordinate[2]],
 [max_coordinate[0], max_coordinate[1], min_coordinate[2]],
 ],
 coords={"corner": range(8), "axis": list(axes)},
 corners = np.transpose(corners, (2, 0, 1))
 else:
 # 3D bounding cube
 assert len(axes) == 3
 corners = np.array(
 [
 [min_coordinate[:, 0], min_coordinate[:, 1], min_coordinate[:, 2]],
 [min_coordinate[:, 0], min_coordinate[:, 1], max_coordinate[:, 2]],
 [min_coordinate[:, 0], max_coordinate[:, 1], max_coordinate[:, 2]],
 [min_coordinate[:, 0], max_coordinate[:, 1], min_coordinate[:, 2]],
 [max_coordinate[:, 0], min_coordinate[:, 1], min_coordinate[:, 2]],
 [max_coordinate[:, 0], min_coordinate[:, 1], max_coordinate[:, 2]],
 [max_coordinate[:, 0], max_coordinate[:, 1], max_coordinate[:, 2]],
 [max_coordinate[:, 0], max_coordinate[:, 1], min_coordinate[:, 2]],
 ]
 )
 corners = np.transpose(corners, (2, 0, 1))
 output = DataArray(
 corners,
 coords={
 "box": range(num_boxes),
 "corner": range(corners.shape[1]),
 "axis": list(axes),
 },
 )
 if num_boxes > 1:
 return output
 return output.squeeze().drop_vars("box")


@nb.njit(parallel=False, nopython=True)
def _create_slices_and_translation(
 min_values: nb.types.Array,
 max_values: nb.types.Array,
) -> tuple[nb.types.Array, nb.types.Array]:
 n_boxes, n_dims = min_values.shape
 slices = np.empty((n_boxes, n_dims, 2), dtype=np.float64) # (n_boxes, n_dims, [min, max])
 translation_vectors = np.empty((n_boxes, n_dims), dtype=np.float64) # (n_boxes, n_dims)

 for i in range(n_boxes):
 for j in range(n_dims):
 slices[i, j, 0] = min_values[i, j]
 slices[i, j, 1] = max_values[i, j]
 translation_vectors[i, j] = np.ceil(max(min_values[i, j], 0))

 return slices, translation_vectors


def _process_data_tree_query_result(query_result: DataTree) -> DataTree | None:
 d = {}
 for k, data_tree in query_result.items():
 v = data_tree.values()
 assert len(v) == 1
 xdata = v.__iter__().__next__()
 if 0 in xdata.shape:
 if k == "scale0":
 return None
 else:
 d[k] = xdata

 # Remove scales after finding a missing scale
 scales_to_keep = []
 for i, scale_name in enumerate(d.keys()):
 if scale_name == f"scale{i}":
 scales_to_keep.append(scale_name)
 else:
 break

 # Case in which scale0 is not present but other scales are
 if len(scales_to_keep) == 0:
 return None

 d = {k: d[k] for k in scales_to_keep}
 result = DataTree.from_dict(d)

 # Rechunk the data to avoid irregular chunks
 for scale in result:
 result[scale]["image"] = result[scale]["image"].chunk("auto")

 return result


def _process_query_result(
 result: DataArray | DataTree, translation_vector: ArrayLike, axes: tuple[str, ...]
) -> DataArray | DataTree | None:
 from spatialdata.transformations import get_transformation, set_transformation

 if isinstance(result, DataArray):
 if 0 in result.shape:
 return None
 # rechunk the data to avoid irregular chunks
 result = result.chunk("auto")
 elif isinstance(result, DataTree):
 result = _process_data_tree_query_result(result)
 if result is None:
 return None

 result = compute_coordinates(result)

 if not np.allclose(np.array(translation_vector), 0):
 translation_transform = Translation(translation=translation_vector, axes=axes)

 transformations = get_transformation(result, get_all=True)
 assert isinstance(transformations, dict)

 new_transformations = {}
 for coordinate_system, initial_transform in transformations.items():
 new_transformation: BaseTransformation = Sequence(
 [translation_transform, initial_transform],
 )
 new_transformations[coordinate_system] = new_transformation
 set_transformation(result, new_transformations, set_all=True)

 # let's make a copy of the transformations so that we don't modify the original object
 t = get_transformation(result, get_all=True)
 assert isinstance(t, dict)
 set_transformation(result, t.copy(), set_all=True)

 return result


def _get_filtered_or_unfiltered_tables(
Expand Down
Loading