File wellmet/voronoi.py changed (mode: 100644) (index 470aebb..4554c99) |
... |
... |
import scipy.stats as stats |
8 |
8 |
from scipy.spatial import KDTree |
from scipy.spatial import KDTree |
9 |
9 |
from scipy.spatial import Delaunay |
from scipy.spatial import Delaunay |
10 |
10 |
from scipy import spatial |
from scipy import spatial |
|
11 |
|
from scipy.spatial.distance import cdist |
11 |
12 |
from scipy import interpolate |
from scipy import interpolate |
12 |
13 |
|
|
13 |
14 |
import collections # for defaultdict, namedtuple |
import collections # for defaultdict, namedtuple |
|
... |
... |
class TwoHorses: |
513 |
514 |
|
|
514 |
515 |
|
|
515 |
516 |
|
|
|
517 |
|
class TwoPoints: |
|
518 |
|
def __init__(self, points, couple_indices): |
|
519 |
|
i, j = couple_indices |
|
520 |
|
|
|
521 |
|
self.point_i = point_i = points[i] |
|
522 |
|
self.point_j = point_j = points[j] |
|
523 |
|
#č vektor zadáme přímkou mezí vzorky |
|
524 |
|
#č od "i" k "j", normalizovat nebudeme |
|
525 |
|
self.normal_vector = point_j - point_i |
|
526 |
|
|
|
527 |
|
#č jako indice zkusme použit bod na usečce uprostřed mezí vzorky |
|
528 |
|
self.half_point = np.mean(points[[i, j]], axis=0) |
|
529 |
|
|
|
530 |
|
self.delimit = np.inner(self.half_point, self.normal_vector) |
|
531 |
|
|
|
532 |
|
#č pokud skalarní součin je menší jak delimit |
|
533 |
|
#č (protože směr normály od "i" k "j"), |
|
534 |
|
#č tak patří hledané do poloprostoru "i" |
|
535 |
|
mask = (points @ self.normal_vector) < self.delimit |
|
536 |
|
|
|
537 |
|
|
|
538 |
|
self.tree_i = KDTree(points[mask]) |
|
539 |
|
self.tree_j = KDTree(points[~mask]) |
|
540 |
|
#č neznám convinient cestu pro získaní |
|
541 |
|
#č tamtamtoho indexu na maskovaném poli |
|
542 |
|
#č np.ma mně příjde zbytečným |
|
543 |
|
# ha! monkey patching |
|
544 |
|
self.tree_i.idx = len(mask[:i][mask[:i]]) |
|
545 |
|
self.tree_j.idx = len(mask[:j][~mask[:j]]) |
|
546 |
|
|
|
547 |
|
#č bacha na indexy. Musíme vyhodit i-tý a j-tý vzorky |
|
548 |
|
mask[i] = False |
|
549 |
|
self.points_i = points[mask] |
|
550 |
|
#č teď ale jedeme se stejnou, ale invertovanou maskou |
|
551 |
|
mask[i] = True |
|
552 |
|
mask[j] = True |
|
553 |
|
self.points_j = points[~mask] |
|
554 |
|
|
|
555 |
|
|
|
556 |
|
def query(self, nodes, order='F'): |
|
557 |
|
#č pokud skalarní součin je menší jak delimit |
|
558 |
|
#č (protože směr normály od "i" k "j"), |
|
559 |
|
#č tak patří hledané do poloprostoru "i" |
|
560 |
|
halfmask = (nodes @ self.normal_vector) < self.delimit |
|
561 |
|
|
|
562 |
|
nodes_i = nodes[halfmask] |
|
563 |
|
nodes_j = nodes[~halfmask] |
|
564 |
|
|
|
565 |
|
|
|
566 |
|
status_i = status_j = False |
|
567 |
|
if nodes_i.size: |
|
568 |
|
status_i, mask_i, dd_i = self._run(self.point_j, self.points_i, |
|
569 |
|
self.tree_j, nodes_i, order) |
|
570 |
|
|
|
571 |
|
if nodes_j.size: |
|
572 |
|
status_j, mask_j, dd_j = self._run(self.point_i, self.points_j, |
|
573 |
|
self.tree_i, nodes_j, order) |
|
574 |
|
|
|
575 |
|
if not (status_i or status_j): |
|
576 |
|
return False, None, None |
|
577 |
|
|
|
578 |
|
#č globální maska |
|
579 |
|
gmask = np.empty(len(nodes), dtype=bool) |
|
580 |
|
|
|
581 |
|
if status_i and not status_j: |
|
582 |
|
gmask[halfmask] = mask_i |
|
583 |
|
gmask[~halfmask] = False |
|
584 |
|
return True, gmask, TData(X=nodes[gmask], imask=mask_i[mask_i], dd=dd_i) |
|
585 |
|
|
|
586 |
|
if status_j and not status_i: |
|
587 |
|
gmask[halfmask] = False |
|
588 |
|
gmask[~halfmask] = mask_j |
|
589 |
|
return True, gmask, TData(X=nodes[gmask], imask=~mask_j[mask_j], dd=dd_j) |
|
590 |
|
|
|
591 |
|
|
|
592 |
|
gmask[halfmask] = mask_i |
|
593 |
|
gmask[~halfmask] = mask_j |
|
594 |
|
|
|
595 |
|
dd_mask = halfmask[gmask] |
|
596 |
|
dd = np.empty(len(dd_i) + len(dd_j)) |
|
597 |
|
dd[dd_mask] = dd_i |
|
598 |
|
dd[~dd_mask] = dd_j |
|
599 |
|
|
|
600 |
|
return True, gmask, TData(X=nodes[gmask], imask=dd_mask, dd=dd) |
|
601 |
|
|
|
602 |
|
|
|
603 |
|
@staticmethod |
|
604 |
|
def _run(center_2, points_1, tree_2, nodes_1, order): |
|
605 |
|
""" |
|
606 |
|
č Úkolem první fáze je vyhodit co nejvíc bodíků |
|
607 |
|
č Bodíky bereme z i-tého poloprostoru a kontrolní vzdálenosti - |
|
608 |
|
č vůči vzorku "j" z protílehlé strany. |
|
609 |
|
č Předpokládáme, že na konci první fázi získame už |
|
610 |
|
č dobře vyfiltrované tečičky, které přece pro formu |
|
611 |
|
č musíme zkontrolovat i vůči vzorkům na opačné straně. |
|
612 |
|
č Předpokladáme, že zůstaly jíž vícemené dobré bodíky |
|
613 |
|
č a že smyčka druhé fáze už by musela běžet přes skoro všechno. |
|
614 |
|
č Proto se nebudeme snažit dohnat Cečkovej kód v Pythonu |
|
615 |
|
č a poslední část pokorně svěříme KDTree - |
|
616 |
|
č nejrychlejší implementaci the nearest search na dívokém západě. |
|
617 |
|
""" |
|
618 |
|
#č cdist vysloveně chce 2D pola, viděl jsem to i v kódě. |
|
619 |
|
#č Vrací sloupcovou matici. Co kdyby s F pořadí by to jelo rychlejc? |
|
620 |
|
#d1 = np.empty((len(nodes), 1), dtype=float, order=order) |
|
621 |
|
#cdist(nodes, [center_1], 'sqeuclidean', out=d1) |
|
622 |
|
d2 = np.empty((len(nodes_1), 1), dtype=float, order=order) |
|
623 |
|
cdist(nodes_1, [center_2], 'sqeuclidean', out=d2) |
|
624 |
|
|
|
625 |
|
mask = np.isfinite(d2) |
|
626 |
|
|
|
627 |
|
# #č udělám to takhle i když neumím si představit, |
|
628 |
|
# #č jak se můžou objevit NaNs a nekonečna ve d2, když nejsou v d1 |
|
629 |
|
# mask = np.isfinite(d1) & np.isfinite(d2) |
|
630 |
|
|
|
631 |
|
n = len(mask[mask]) |
|
632 |
|
d_i = np.empty((n, 1), dtype=float, order=order) |
|
633 |
|
|
|
634 |
|
#č pro nás spíše budou kritické poslední přídané vzorky |
|
635 |
|
for point in reversed(points_1): |
|
636 |
|
cdist(nodes_1[mask[:,0]], np.atleast_2d(point), 'sqeuclidean', out=d_i) |
|
637 |
|
|
|
638 |
|
mask[mask[:,0]] &= d2[mask[:,0]] <= d_i |
|
639 |
|
|
|
640 |
|
n = len(mask[mask]) |
|
641 |
|
if n == 0: |
|
642 |
|
return False, None, None #mask[:,0], d_2[mask_2] |
|
643 |
|
d_i.resize((n, 1), refcheck=False) |
|
644 |
|
|
|
645 |
|
#č ty špinavé 2D matice zhora jíž nepotřebujem |
|
646 |
|
#č převedeme (co nejlevnějc) masku na 1D |
|
647 |
|
mask = mask[:,0] |
|
648 |
|
d_2, j = tree_2.query(nodes_1[mask], k=1, p=2, workers=-1) |
|
649 |
|
mask_2 = (j == tree_2.idx) |
|
650 |
|
|
|
651 |
|
if np.any(mask_2): |
|
652 |
|
#č hotovo. Vzdálenosti jíž nemusíme kontrolovat. |
|
653 |
|
#č ve smyčce d2 jsme jíž prohnali. |
|
654 |
|
mask[mask] = mask_2 |
|
655 |
|
return True, mask, d_2[mask_2] |
|
656 |
|
else: |
|
657 |
|
return False, None, None |
|
658 |
|
|
|
659 |
|
|
|
660 |
|
|
|
661 |
|
# @staticmethod |
|
662 |
|
# def _first_query(center, points, nodes, order): |
|
663 |
|
# |
|
664 |
|
# d1 = np.empty((len(nodes), 1), dtype=float, order=order) |
|
665 |
|
# d2 = np.full((len(nodes), 1), np.inf, dtype=float, order=order) |
|
666 |
|
# cdist(nodes, [center], 'sqeuclidean', out=d1) |
|
667 |
|
# mask = np.isfinite(d1) |
|
668 |
|
# |
|
669 |
|
# n = len(mask[mask]) |
|
670 |
|
# #d2 = np.full((n, 1), np.inf, dtype=float, order=order) |
|
671 |
|
# d3 = np.empty((n, 1), dtype=float, order=order) |
|
672 |
|
# |
|
673 |
|
# for point in reversed(points): |
|
674 |
|
# cdist(nodes[mask[:,0]], np.atleast_2d(point), 'sqeuclidean', out=d3) |
|
675 |
|
# |
|
676 |
|
# np.fmin(d2[mask[:,0]], d3, out=d2[mask[:,0]]) |
|
677 |
|
# mask[mask[:,0]] &= d1[mask[:,0]] <= d3 |
|
678 |
|
# |
|
679 |
|
# n = len(mask[mask]) |
|
680 |
|
# if n == 0: |
|
681 |
|
# break |
|
682 |
|
# d3.resize((n, 1), refcheck=False) |
|
683 |
|
# return mask, d1, d2 |
|
684 |
|
# |
|
685 |
|
# @staticmethod |
|
686 |
|
# def _second_query(center, points, nodes, order): |
|
687 |
|
# |
|
688 |
|
# d1 = np.empty((len(nodes), 1), dtype=float, order=order) |
|
689 |
|
# cdist(nodes, [center], 'sqeuclidean', out=d1) |
|
690 |
|
# mask = np.isfinite(d1) |
|
691 |
|
# n = len(mask[mask]) |
|
692 |
|
# d3 = np.empty((n, 1), dtype=float, order=order) |
|
693 |
|
# |
|
694 |
|
# for point in reversed(points): |
|
695 |
|
# cdist(nodes[mask[:,0]], np.atleast_2d(point), 'sqeuclidean', out=d3) |
|
696 |
|
# |
|
697 |
|
# mask[mask[:,0]] &= d1[mask[:,0]] <= d3 |
|
698 |
|
# |
|
699 |
|
# n = len(mask[mask]) |
|
700 |
|
# if n == 0: |
|
701 |
|
# break |
|
702 |
|
# d3.resize((n, 1), refcheck=False) |
|
703 |
|
# return mask, d1 |
|
704 |
|
|
|
705 |
|
|
|
706 |
|
|
|
707 |
|
|
|
708 |
|
|
|
709 |
|
|
516 |
710 |
|
|
517 |
711 |
#оӵ кык точкаен Вороной |
#оӵ кык точкаен Вороной |
518 |
712 |
# KickTouchCayenneVoronoi |
# KickTouchCayenneVoronoi |