Commit 98738e22 authored by casagran's avatar casagran
Browse files

backup stuff

parent cf0690e0
/home/experience/psiche/com-psiche/2021/Proudhon_0621
\ No newline at end of file
......@@ -833,7 +833,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.4"
"version": "3.9.5"
}
},
"nbformat": 4,
......
1st step: go to bigmeca
option 1: port forward (prefer this if the second is too slow)
ssh -L 8888:localhost:8888 com-psiche@195.221.7.39
don't use -X because it might conflict with someone else
not critical
"-L (localhost:8888):localhost:8888"
tells ssh to bridge the port 8888 in the local machine with the port 8888 of the server
195.221.7.39
bigmeca's ip address
com-psiche
user ==> necessary to access the data in /data/experiences/psiche/com-psiche
option 2: ssh with graphics (using a firefox in the server machine and sending the graphics)
ssh -C -X com-psiche@195.221.7.39
-C s
###If not need ssh start here
2nd step: launch jupyter lab
if option1: open a browser on your local machine
if option2: open a browser on the server machine using the ssh-connected terminal
it's better if you create a tmux so that other stuff don't eventually break
to debug if jupyter lab doesnt open, check that it is running in tmux
check open session:
tmux list-session
open the one for jupyter lab
tmux attach-session -t jlab
"-t" for "target"
to quit: Ctrl+b, d
'Ctrl+b' is the prefix to any command with tmux
"d" for detach
or create a new one
tmux new-session -n jlab
"-n" for name
if you indeed need to relaunch jupyter lab, use the correct conda env
see the available conda envs
conda env list
there should be a 'jlab' or 'jupyterlab'
the current env has a '*'
activate it
conda activate jlab
then launch launch jupyter lab
jupyter lab
quit your tmux (Ctrl+b, d)
recheck your local browser
where pydct is
remember that dct is a repo
./dct/zUtil_Python/pydct/
full: /home/experiences/psiche/com-psiche/dct/zUtil_Python/pydct/
env you're supposed to use: pydct-dev
the envs are in: pydct/envs
\ No newline at end of file
import logging
from pathlib import Path
import matplotlib
from IPython.core.interactiveshell import InteractiveShell
from pydct import common, log, parse
_logger = logging.getLogger("pydct.segment.config")
_logger.handlers = []
_logger.propagate = False
_logger.addHandler(log.get_stdout_handler())
_logger.setLevel(logging.INFO)
this_dir = Path(__file__).parent
pydct_src_dir = this_dir.parent
matplotlib.style.use(pydct_src_dir / "pydct.mplstyle")
InteractiveShell.ast_node_interactivity = "all"
import functools
import logging
from enum import Enum
import ipywidgets as widgets
import matplotlib.pyplot as plt
import numpy as np
import skimage
from IPython import display
from ipywidgets import fixed, interact, interact_manual, interactive
from matplotlib import pyplot as plt
from numpy import ndarray
from pydct import common, log, parse
from skimage import filters, measure
_logger = logging.getLogger("pydct.segment")
_logger.handlers = []
_logger.propagate = False
_logger.addHandler(log.get_stdout_handler())
_logger.setLevel(logging.INFO)
# threshold = .6
# min_3dblob_volume = 10
def generate_seeds_binary_threshold(
data: ndarray, threshold: float, min_3dblob_volume: int, return_labels_volume=False
) -> (ndarray, int):
_logger.debug(f"applying threshold \\approx {float(f'{threshold:.4g}'):,.3f}")
# isfinite is necessary because we exclude the direct beam with nan
seeds = (data >= threshold) & np.isfinite(data)
_logger.debug(f"finding connected components")
# each connected region gets its own label (integer)
labels = skimage.measure.label(seeds)
_logger.debug(f"getting region props")
# a list of volumes of the blobs in the same order as their labels
props = skimage.measure.regionprops(labels, data, cache=False)
if not return_labels_volume:
_logger.debug("getting rid of the labels volume")
del labels
volumes = np.array([p.area for p in props])
_logger.debug(f"{len(props)} found, now filtering them by volume")
# filter the labels of regions with enough volume
selected_labels = volumes > min_3dblob_volume
_logger.debug(f"{np.sum(selected_labels)} left, now making getting a list of the label values left")
selected_labels_values = np.array([p.label for p in props])[selected_labels]
if return_labels_volume:
_logger.debug(f"erasing the discarde labeled regions from the labels volume")
# mask the labels volume filtering those without enough volume
selected_seeds = labels * np.isin(labels, selected_labels_values).astype(np.int32)
_logger.debug(f"now filtering the region props")
selected_props = [p for p in props if p.label in selected_labels_values]
if return_labels_volume:
return selected_seeds, selected_props, selected_labels.size
else:
return None, selected_props, selected_labels.size
from typing import List
from pandas import DataFrame
from skimage.measure._regionprops import RegionProperties
def get_seed_region_max(data: ndarray, rp: RegionProperties):
# find max value in the seed blob, then its absolute position and value
# 'absolute position' = in the whole image's frame
# 'relative position' = in the subvolume ("ROI") around the seed's maximum
# at this moment the ROI == BB (bounding box) of the blob
# but later it will not necessarily be
rel_position = np.unravel_index(
np.nanargmax(rp.intensity_image), rp.intensity_image.shape # serial index, not multi-dimensioned
)
# origin in the blob's referential
abs_bb_origin = rp.bbox[0:3]
abs_max_position = np.array(abs_bb_origin) + np.array(rel_position)
max_value = data[tuple(abs_max_position)]
return dict(
x=abs_max_position[1],
y=abs_max_position[2],
z=abs_max_position[0],
v=max_value,
# these will serve to filter the seeds as we find bad ones
bad=False,
label=rp.label,
)
from dataclasses import dataclass
from typing import Tuple
@dataclass
class DiffSpot:
# the bounding box in the image's indexing system
# using the 3d upper left and 3d lower right corners
bb: Tuple[Tuple[int, int, int], Tuple[int, int, int]]
mask: ndarray
def _worker_init(
shared_raw_array,
data_shape,
data_dtype,
):
global global_data_
global_data_ = np.frombuffer(shared_raw_array, dtype=data_dtype).reshape(data_shape)
from enum import Enum
class RejectionCode(Enum):
too_small = 0
too_big = 1
seed_max_is_not_max = 2
from dataclasses import dataclass
@dataclass
class Rejection:
code: RejectionCode
info: dict
def _grow_one_seed_no_iteration(
abs_max_position,
max_value,
data_,
threshold_ratio_,
blob_max_shape_,
blob_min_shape_,
is_parallel=False,
return_rejections=False,
):
if is_parallel:
global global_data_
data_ = global_data_
assert data_ is not None
# ~local~ threshold (based on the seed's maximum)
threshold = threshold_ratio_ * max_value
# the bounding box here is defined by the upper left (ul)
# and bottom right (br) corners,
# where the ul (ul.z, ul.x, ul.y) is inclusive but not the br (br.z, br.x, br.y)
# so one can range(ul[0], br[0])
# as in skimage
# 3D: i use 3 axes here but the naming is based on the 2D case ("upper left"...)
max_roi = np.stack(
[abs_max_position - (blob_max_shape_ // 2), abs_max_position + (blob_max_shape_ // 2) + (blob_max_shape_ % 2)]
)
# make sure the UL/BR corners of the ROI don't go out of the image
max_roi = np.array(
[
# ul
[max(max_roi[0, ax], 0) for ax in range(3)],
# br
[min(max_roi[1, ax], data_.shape[ax]) for ax in range(3)],
]
)
# we first assume the max's position to be in the center of the roi
# but if it's close to the image's borders, the effective roi's
# size might get smaller due to the previous operation, so we
# compensate it by extending it again
for ax in range(3):
# roi touches the image's lower bound
if max_roi[0, ax] == 0:
# throw the opposite corner to the maximum
max_roi[1, ax] = min(blob_max_shape_[ax], data_.shape[ax])
# roi touches the image's upper bound
elif max_roi[1, ax] == data_.shape[ax]:
# do the opposite
max_roi[0, ax] = max(data_.shape[ax] - blob_max_shape_[ax], 0)
seed_max_rel_position = abs_max_position - max_roi[0]
# cutoff the roi
slicez = slice(*max_roi[:, 0])
slicex = slice(*max_roi[:, 1])
slicey = slice(*max_roi[:, 2])
roi_data = data_[slicez, slicex, slicey]
# find the seed's CC region
roi_labels = skimage.measure.label(roi_data >= threshold)
seed_max_label = roi_labels[tuple(seed_max_rel_position)]
roi_mask = roi_labels == seed_max_label
roi_data_masked = (roi_data * roi_mask.astype(np.float32)).astype(data_.dtype)
max_is_from_seed = max_value >= np.nanmax(roi_data_masked)
if not max_is_from_seed:
if not return_rejections:
return
return Rejection(
code=RejectionCode.seed_max_is_not_max,
info=dict(
abs_max_position=abs_max_position,
max_value=max_value,
threshold=threshold,
argmax=np.nanargmax(roi_data_masked),
),
)
# there should only be one because the mask got rid of others,
# and at least one because the process starts with one
grown_blob_region_props = skimage.measure.regionprops(roi_mask.astype(np.int32), cache=False)[0]
# bounding box relative to the ROI
grown_blob_relative_bb = np.array(grown_blob_region_props.bbox).reshape((2, 3))
grown_blob_bb_shape = np.array(grown_blob_region_props.image.shape)
# if the bbox relative to the roi is the same size as the roi in any dimension
# then it is touching both borders ==> too big
is_too_big = np.any(grown_blob_bb_shape >= blob_max_shape_)
# ul = upper left corner
# br = bottom right corner
ul = grown_blob_relative_bb[0] + max_roi[0]
br = grown_blob_relative_bb[1] + max_roi[0]
abs_bounding_box = (tuple(ul), tuple(br))
if is_too_big:
if not return_rejections:
return
return Rejection(
code=RejectionCode.too_big,
info=dict(
abs_max_position=abs_max_position,
max_value=max_value,
threshold=threshold,
grown_blob_bb_shape=grown_blob_bb_shape,
grown_blob_region_props=grown_blob_region_props,
abs_bounding_box=abs_bounding_box,
),
)
is_too_small = np.any(grown_blob_bb_shape < blob_min_shape_)
if is_too_small:
if not return_rejections:
return
return Rejection(
code=RejectionCode.too_small,
info=dict(
abs_max_position=abs_max_position,
max_value=max_value,
threshold=threshold,
grown_blob_bb_shape=grown_blob_bb_shape,
grown_blob_region_props=grown_blob_region_props,
abs_bounding_box=abs_bounding_box,
),
)
return DiffSpot(bb=abs_bounding_box, mask=grown_blob_region_props.image.astype(np.int8))
import contextlib
from multiprocessing import Pool
def grow_seeds_no_iteration_parallel(
data_,
shared_raw_array_,
seed_maxima_: list,
threshold_ratio_,
blob_min_shape_,
blob_max_shape_,
nprocs=1,
return_rejections=False,
):
blob_min_shape_ = np.array(blob_min_shape_)
blob_max_shape_ = np.array(blob_max_shape_)
maxima_df = DataFrame.from_records(data=seed_maxima_).sort_values(by=["v"], ascending=False)
blobs = list()
rejections = dict()
if nprocs == 1:
for idx in range(maxima_df.shape[0]):
row = maxima_df.iloc[idx]
# row.v = the maximum value of the blob
value = row.v
position = row[["z", "x", "y"]].values.astype(np.int32)
blob = _grow_one_seed_no_iteration(
position,
value,
data_,
threshold_ratio_,
blob_max_shape_,
blob_min_shape_,
is_parallel=False,
return_rejections=return_rejections,
)
if isinstance(blob, DiffSpot):
blobs.append(blob)
elif return_rejections:
blob.info["label"] = row.label
rejections[row.label] = blob
return blobs if not return_rejections else (blobs, rejections)
assert shared_raw_array_ is not None
with contextlib.closing(
Pool(
processes=nprocs,
initializer=_worker_init,
initargs=(
shared_raw_array_,
data_.shape,
data_.dtype,
),
)
) as pool:
blobs = pool.starmap_async(
_grow_one_seed_no_iteration,
(
(
maxima_df.iloc[idx][["z", "x", "y"]].values.astype(np.int32), # position
maxima_df.iloc[idx]["v"], # value
None, # data
threshold_ratio_,
blob_max_shape_,
blob_min_shape_,
True, # is_parallel
False, # return_rejections
)
for idx in range(maxima_df.shape[0])
),
)
pool.join()
assert blobs.ready() and blobs.successful()
return [b for b in blobs.get() if b is not None]
import sys
def memory_usage(dir_return_value):
"""
src: https://stackoverflow.com/questions/1823058/how-to-print-number-with-commas-as-thousands-separators
"""
global globals
# These are the usual ipython objects, including this one you are creating
ipython_vars = ["In", "Out", "exit", "quit", "get_ipython", "ipython_vars"]
# Get a sorted list of the objects and their sizes
return sorted(
[
(x, id(globals().get(x)), f"{sys.getsizeof(globals().get(x)) / 1024**2:,.3f} Mb")
for x in dir_return_value
if not x.startswith("_") and x not in sys.modules and x not in ipython_vars
],
key=lambda x: x[2],
reverse=True,
)
from typing import List
import ipywidgets
import numpy as np
from IPython.display import display
from ipywidgets import HBox, IntSlider, VBox, interact, link
from matplotlib import pyplot as plt
from numpy import ndarray
def imshow3d(img3d: ndarray, zidx: int, ax=None, **imshow_kwargs):
assert img3d.ndim == 3
nz, nx, ny = img3d.shape
assert 0 <= zidx < nz
if ax is None:
h = 5 # height
w = h * nx / ny # width
fig, ax = plt.subplots(1, 1, figsize=(h, w))
# this has to be set like this so that the colors
# of different frames will be consistent
# otherwise the min/max are reset at each frame
try:
vmin = imshow_kwargs["vmin"]
except KeyError:
vmin = np.nanmin(img3d)
try:
vmax = imshow_kwargs["vmax"]
except KeyError:
vmax = np.nanmax(img3d)
# replace defaults with user-given values
imshow_kwargs = {
**dict(
vmin=vmin,
vmax=vmax,
cmap="gray",
),
**imshow_kwargs,
}
# just a hot fix...
if "norm" in imshow_kwargs:
del imshow_kwargs['vmin']
del imshow_kwargs['vmax']
im = ax.imshow(img3d[zidx, :, :], **imshow_kwargs)
return im
def get_imshow3d_widget(data3d, ax=None, **imshow_kwargs):
assert data3d.ndim == 3
nz, nx, ny = data3d.shape