Valis blocks at "Measuring error"

[marco1108]

Hi all,

I am running Valis on my MacBook Pro in a conda environment. I do not get any error message, it just runs for hours the measuring error until I terminate the script.

Has anyone encountered this before?

Thanks

M

python /Users/marcopisa/Library/CloudStorage/OneDrive-Nexus365/Projects/BAPP/Valis_Script/VALIS_modified.py $DIRECTORY $PREFIX

done
1
/Users/marcopisa/Library/CloudStorage/OneDrive-Nexus365/Projects/BAPP/Valis/MS363C
MS363_SC_C
['/Users/marcopisa/Library/CloudStorage/OneDrive-Nexus365/Projects/BAPP/Valis_Script/VALIS_modified.py', '/Users/marcopisa/Library/CloudStorage/OneDrive-Nexus365/Projects/BAPP/Valis/MS363C', 'MS363_SC_C']
MS363_SC_C_fibrinogen_40X.svs
SLF4J: Failed to load class "org.slf4j.impl.StaticLoggerBinder".
SLF4J: Defaulting to no-operation (NOP) logger implementation
SLF4J: See http://www.slf4j.org/codes.html#StaticLoggerBinder for further details.
JVM has been initialized. Be sure to call registration.kill_jvm() or slide_io.kill_jvm() at the end of your script.
/Users/marcopisa/opt/anaconda3/envs/valis-env/lib/python3.10/site-packages/valis/valtils.py:21: UserWarning: found 3 matches for MS363C: /Users/marcopisa/Library/CloudStorage/OneDrive-Nexus365/Projects/BAPP/Valis/MS363C/MS363_SC_C_BAPP_40X.svs, /Users/marcopisa/Library/CloudStorage/OneDrive-Nexus365/Projects/BAPP/Valis/MS363C/MS363_SC_C_40X.geojson, /Users/marcopisa/Library/CloudStorage/OneDrive-Nexus365/Projects/BAPP/Valis/MS363C/MS363_SC_C_fibrinogen_40X.svs
warnings.warn(warning_msg, warning_type)

==== Converting images

100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 2/2 [00:03<00:00, 1.56s/it]

==== Processing images

100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 2/2 [00:03<00:00, 1.90s/it]

==== Normalizing images

100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 2/2 [00:00<00:00, 20.30it/s]

==== Rigid registraration

100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 2/2 [00:00<00:00, 8.81it/s]

======== Detecting features

100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 2/2 [00:00<00:00, 4.85it/s]

======== Matching images

100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 1.49it/s]

======== Sorting images

======== Calculating transformations

100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 2605.16it/s]
100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 2/2 [00:00<00:00, 7436.71it/s]

======== Rigid registration complete in 1.199 seconds

==== Non-rigid registraration

======== Preparing images for non-rigid registration

100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 2/2 [00:01<00:00, 1.69it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 2/2 [00:04<00:00, 2.29s/it]

======== Non-rigid registration complete in 5.276 seconds

==== Measuring error

100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 2/2 [00:00<00:00, 2.92it/s]
2305it [08:24, 4.59it/s]

[cdgatenbee]

Hi @marco1108,
Thanks for reporting this. It is strange that valis would get hung up measuring the error, since it's a just few basic calculations. Could you let me know if results are getting saved to the data directory? If so, that would indicate that Valis.measure_error has completed, and that valis is getting hung up somewhere after that. Could you confirm that you are shutting down the JVM with registration.kill_jvm() or slide_io.kill_jvm()? If not, that could be what is preventing the run from completing, since Python won't stop running until the JVM gets shutdown. If you are shutting down the JVM, could you please share the code you are using to run valis?

Best,
-Chandler

[marco1108]

Hi @cdgatenbee,

Thanks for you prompt reply.

I managed to understand why the script was blocked in "measuring error" stage. The .geojson file I was using contained not only the information on the annotations (that I wish to transpose to consecutive tissue sections thorugh Valis) but also all the positive cell detections. This is why Valis was unable to complete task even after many hours, it was trying to transpose all the cell detections.

I have now deleted all the cell detections from the .geojson file, and the script works without errors. The script has generated a new folder that contains all the co-registratio files. From review of these outputs it seems that Valis has done an amazing job in co-registering the consecutive images.

However, I cannot find the output that I needed which is a new .geogeson with the adapted annotations, so that I can use them on the consecutive slides Valis has co-registered (or alternatively, a new .svs file with the co-registered image - as I now only have a low resolution .png).

This is the new folder created by the script:

Here is the script used:
`import numpy as np
import json
import os
import pathlib
import shapely
import tqdm
from copy import deepcopy

from skimage import exposure
import numpy as np
from valis import preprocessing

class PreProcessor(preprocessing.ImageProcesser):
"""Preprocess images for registration

Uses distance to background color and also subtracts the background intensity

"""

def __init__(self, image, src_f, level, series, *args, **kwargs):
    super().__init__(image=image, src_f=src_f, level=level,
                     series=series, *args, **kwargs)


def create_mask(self):
    _, tissue_mask = preprocessing.create_tissue_mask_from_rgb(self.image)

    return tissue_mask


def process_image(self, brightness_q=0.99, *args, **kwargs):

    processed_img, cam = preprocessing.calc_background_color_dist(self.image, brightness_q=brightness_q)
    processed_img = exposure.rescale_intensity(processed_img, in_range="image", out_range=(0, 1))

    # Subtract background
    brightest_thresh = np.quantile(cam[..., 0], 0.9)
    brightest_idx = np.where(cam[..., 0] >= brightest_thresh)
    brightest_pixels = processed_img[brightest_idx]
    brightest_rgb = np.median(brightest_pixels, axis=0)
    no_bg = processed_img - brightest_rgb
    no_bg = np.clip(no_bg, 0, 255)

    # Adjust range and perform adaptive histogram equalization
    no_bg = (255*exposure.equalize_adapthist(no_bg/no_bg.max())).astype(np.uint8)

    return no_bg

from valis import registration, valtils, warp_tools

def clip_xy(xy, shape_rc):
"""Clip xy coordintaes to be within image

"""
clipped_x =  np.clip(xy[:, 0], 0, shape_rc[1])
clipped_y =  np.clip(xy[:, 1], 0, shape_rc[0])

clipped_xy = np.dstack([clipped_x, clipped_y])[0]
return clipped_xy

def _warp_shapely(geom, warp_fxn, warp_kwargs):
"""Warp a shapely geometry
Based on shapely.ops.trasform

"""
if "dst_shape_rc" in warp_kwargs:
    dst_shape_rc = warp_kwargs["dst_shape_rc"]
elif "to_dst_shape_rc" in warp_kwargs:
    dst_shape_rc = warp_kwargs["to_dst_shape_rc"]
else:
    dst_shape_rc  = None


if geom.geom_type in ("Point", "LineString", "LinearRing", "Polygon"):
    if geom.geom_type in ("Point", "LineString", "LinearRing"):
        warped_xy = warp_fxn(np.vstack(geom.coords), **warp_kwargs)
        if dst_shape_rc is not None:
            warped_xy = clip_xy(warped_xy, dst_shape_rc)

        return type(geom)(warped_xy.tolist())

    elif geom.geom_type == "Polygon":
        shell_xy = warp_fxn(np.vstack(geom.exterior.coords), **warp_kwargs)
        if dst_shape_rc is not None:
            shell_xy = clip_xy(shell_xy, dst_shape_rc)
        shell = type(geom.exterior)(shell_xy.tolist())
        holes = []
        for ring in geom.interiors:
            holes_xy = warp_fxn(np.vstack(ring.coords), **warp_kwargs)
            if dst_shape_rc is not None:
                holes_xy = clip_xy(holes_xy, dst_shape_rc)

            holes.append(type(ring)(holes_xy))

        return type(geom)(shell, holes)

elif geom.geom_type.startswith("Multi") or geom.geom_type == "GeometryCollection":
    return type(geom)([_warp_shapely(part, warp_fxn, warp_kwargs) for part in geom.geoms])
else:
    raise shapely.errors.GeometryTypeError(f"Type {geom.geom_type!r} not recognized")

def warp_shapely_geom_from_to(geom, from_M=None, from_transformation_src_shape_rc=None,
from_transformation_dst_shape_rc=None, from_src_shape_rc=None,
from_dst_shape_rc=None,from_bk_dxdy=None, from_fwd_dxdy=None,
to_M=None, to_transformation_src_shape_rc=None,
to_transformation_dst_shape_rc=None, to_src_shape_rc=None,
to_dst_shape_rc=None, to_bk_dxdy=None, to_fwd_dxdy=None):
"""
Warp xy points using M and/or bk_dxdy/fwd_dxdy. If bk_dxdy is provided, it will be inverted to create fwd_dxdy

Parameters
----------
geom : shapely.geometery
    Shapely geom to warp

M : ndarray, optional
     3x3 affine transformation matrix to perform rigid warp

transformation_src_shape_rc : (int, int)
    Shape of image that was used to find the transformation.
    For example, this could be the original image in which features were detected

transformation_dst_shape_rc : (int, int), optional
    Shape of the image with shape `transformation_src_shape_rc` after warping.
    This could be the shape of the original image after applying `M`.

src_shape_rc : optional, (int, int)
    Shape of the image from which the points originated. For example,
    this could be a larger/smaller version of the image that was
    used for feature detection.

dst_shape_rc : optional, (int, int)
    Shape of image (with shape `src_shape_rc`) after warping

bk_dxdy : ndarray, pyvips.Image
    (2, N, M) numpy array of pixel displacements in the x and y
    directions from the reference image. dx = bk_dxdy[0],
    and dy=bk_dxdy[1]. If `bk_dxdy` is not None, but
    `fwd_dxdy` is None, then `bk_dxdy` will be inverted to warp `xy`.

fwd_dxdy : ndarray, pyvips.Image
    Inverse of bk_dxdy. dx = fwd_dxdy[0], and dy=fwd_dxdy[1].
    This is what is actually used to warp the points.

pt_buffer : int
    If `bk_dxdy` or `fwd_dxdy` are pyvips.Image object, then
    pt_buffer` determines the size of the window around the point used to
    get the local displacements.


Returns
-------
warped_geom : shapely.geom
   Warped `geom`

"""

warp_kwargs = {"from_M": from_M,
               "from_transformation_src_shape_rc": from_transformation_src_shape_rc,
               "from_transformation_dst_shape_rc": from_transformation_dst_shape_rc,
               "from_src_shape_rc": from_src_shape_rc,
               "from_dst_shape_rc":from_dst_shape_rc,
               "from_bk_dxdy":from_bk_dxdy,
               "from_fwd_dxdy":from_fwd_dxdy,
               "to_M":to_M,
               "to_transformation_src_shape_rc": to_transformation_src_shape_rc,
               "to_transformation_dst_shape_rc": to_transformation_dst_shape_rc,
               "to_src_shape_rc": to_src_shape_rc,
               "to_dst_shape_rc": to_dst_shape_rc, "to_bk_dxdy": to_bk_dxdy,
               "to_fwd_dxdy":to_fwd_dxdy}

warped_geom = _warp_shapely(geom, warp_tools.warp_xy_from_to, warp_kwargs)

return warped_geom

def warp_geojson_from_to(geojson_f, annotation_slide, to_slide_obj, src_slide_level=0, src_pt_level=0,
dst_slide_level=0, non_rigid=True):
"""Warp geoms in geojson file from annotation slide to another unwarped slide

Takes a set of geometries found in this annotation slide, and warps them to
their position in the unwarped "to" slide.

Parameters
----------
geojson_f : str
    Path to geojson file containing the annotation geometries

to_slide_obj : Slide
    Slide to which the points will be warped. I.e. `xy`
    will be warped from this Slide to their position in
    the unwarped slide associated with `to_slide_obj`.

src_pt_level: int, tuple, optional
    Pyramid level of the slide/image in which `xy` originated.
    For example, if `xy` are from the centroids of cell segmentation
    performed on the unwarped full resolution image, this should be 0.
    Alternatively, the value can be a tuple of the image's shape (row, col)
    from which the points came. For example, if `xy` are  bounding
    box coordinates from an analysis on a lower resolution image,
    then pt_level is that lower resolution image's shape (row, col).

dst_slide_level: int, tuple, optional
    Pyramid level of the slide/image in to `xy` will be warped.
    Similar to `src_pt_level`, if `dst_slide_level` is an int then
    the points will be warped to that pyramid level. If `dst_slide_level`
    is the "to" image's shape (row, col), then the points will be warped
    to their location in an image with that same shape.

non_rigid : bool, optional
    Whether or not to conduct non-rigid warping. If False,
    then only a rigid transformation will be applied.

"""

if np.issubdtype(type(src_pt_level), np.integer):
    src_pt_dim_rc = annotation_slide.slide_dimensions_wh[src_pt_level][::-1]
else:
    src_pt_dim_rc = np.array(src_pt_level)

if np.issubdtype(type(dst_slide_level), np.integer):
    to_slide_src_shape_rc = to_slide_obj.slide_dimensions_wh[dst_slide_level][::-1]
else:
    to_slide_src_shape_rc = np.array(dst_slide_level)

if src_slide_level != 0:
    if np.issubdtype(type(src_slide_level), np.integer):
        aligned_slide_shape = annotation_slide.val_obj.get_aligned_slide_shape(src_slide_level)
    else:
        aligned_slide_shape = np.array(src_slide_level)
else:
    aligned_slide_shape = annotation_slide.aligned_slide_shape_rc

if non_rigid:
    src_fwd_dxdy = annotation_slide.fwd_dxdy
    dst_bk_dxdy = to_slide_obj.bk_dxdy

else:
    src_fwd_dxdy = None
    dst_bk_dxdy = None

with open(geojson_f) as f:
    annotation_geojson = json.load(f)

warped_features = [None]*len(annotation_geojson["features"])
for i, ft in tqdm.tqdm(enumerate(annotation_geojson["features"])):
    geom = shapely.geometry.shape(ft["geometry"])
    warped_geom = warp_shapely_geom_from_to(geom=geom,
                                    from_M=annotation_slide.M,
                                    from_transformation_dst_shape_rc=annotation_slide.reg_img_shape_rc,
                                    from_transformation_src_shape_rc=annotation_slide.processed_img_shape_rc,
                                    from_dst_shape_rc=aligned_slide_shape,
                                    from_src_shape_rc=src_pt_dim_rc,
                                    from_fwd_dxdy=src_fwd_dxdy,
                                    to_M=to_slide_obj.M,
                                    to_transformation_src_shape_rc=to_slide_obj.processed_img_shape_rc,
                                    to_transformation_dst_shape_rc=to_slide_obj.reg_img_shape_rc,
                                    to_src_shape_rc=to_slide_src_shape_rc,
                                    to_dst_shape_rc=aligned_slide_shape,
                                    to_bk_dxdy=dst_bk_dxdy
                                    )

    warped_ft = deepcopy(ft)
    warped_ft["geometry"] = shapely.geometry.mapping(warped_geom)
    warped_features[i] = warped_ft

warped_geojson = {"type":annotation_geojson["type"], "features":warped_features}

return warped_geojson

if name == "main":

import sys

print(sys.argv, file=sys.stderr)
if (len(sys.argv) < 3):
    print("""ERROR: 2 arguments needed.

Usage: VALIS_modified.py
Exiting.""", file=sys.stderr)
sys.exit(1)

## FIRST ARG slide, results, annotation dir
slide_src_dir = sys.argv[1]
results_dst_dir = sys.argv[1]
annotation_dir = sys.argv[1]

## SECOND ARG file prefix
prefix = sys.argv[2]

annotation_geojson_f = prefix + "_40X.geojson"
annotation_img_f = prefix + "_BAPP_40X.svs"
target_img_f_fibrinogen = prefix + "_fibrinogen_40X.svs"


## for testing only; remove when done
print(target_img_f_fibrinogen)
## sys.exit(0)
##

# slide_src_dir = "/Users/marcopisa/Library/CloudStorage/OneDrive-Nexus365/Projects/BAPP/Valis/MS371C"
# results_dst_dir = "/Users/marcopisa/Library/CloudStorage/OneDrive-Nexus365/Projects/BAPP/Valis/MS371C"
# annotation_dir = "/Users/marcopisa/Library/CloudStorage/OneDrive-Nexus365/Projects/BAPP/Valis/MS371C"
# annotation_geojson_f = "MS371_SC_C_40X.geojson"
# annotation_img_f = "MS371_SC_C_BAPP_40X.svs"
# target_img_f_fibrinogen = "MS371_SC_C_fibrinogen_40X.svs"


warped_geojson_annotation_f_fibrinogen = f"{valtils.get_name(target_img_f_fibrinogen)}_annotations.geojson"

# Perform registration
registrar = registration.Valis(slide_src_dir, results_dst_dir, reference_img_f=annotation_img_f, align_to_reference=True)
rigid_registrar, non_rigid_registrar, error_df = registrar.register(brightfield_processing_cls=PreProcessor)

# Transfer annotation
annotation_geojson_f = os.path.join(annotation_dir, annotation_geojson_f)
# annotation_geojson_f = os.path.join(annotation_dir, "MS371_SC_C_40X.geojson")
annotation_source_slide = registrar.get_slide(annotation_img_f)
target_slide_fibrinogen = registrar.get_slide(target_img_f_fibrinogen)


warped_geojson_fibrinogen = warp_geojson_from_to(annotation_geojson_f, annotation_source_slide, target_slide_fibrinogen)

# Save annotation, which can be dragged and dropped into QuPath
with open(warped_geojson_annotation_f_fibrinogen, 'w') as f:
    json.dump(warped_geojson_fibrinogen, f)

`

[cdgatenbee]

Hi @marco1108,
Glad to hear you figured out what was preventing the script from completing, and that the registrations are looking good! Based on the code you've shared, it looks like the warped geojson file would be in the same directory as VALIS_modified.py, or wherever it is that you are launching the script. If you'd like warped_geojson_annotation_f_fibrinogen to be saved in the sample folder, you could change

warped_geojson_annotation_f_fibrinogen = f"{valtils.get_name(target_img_f_fibrinogen)}_annotations.geojson"

to

import os

warped_geojson_annotation_f_fibrinogen = os.path.join(results_dst_dir, f{valtils.get_name(target_img_f_fibrinogen)}_annotations.geojson")

which should save the updated results in MS317C/MS317C. Please try that out, and if it doesn't work let me know.

Best,
-Chandler

[marco1108]

Hi Chandler,

Amazing, thanks! It works now!

I still get these error messages at the beginning of the output:
"SLF4J: Failed to load class "org.slf4j.impl.StaticLoggerBinder".
SLF4J: Defaulting to no-operation (NOP) logger implementation
SLF4J: See http://www.slf4j.org/codes.html#StaticLoggerBinder for further details.
JVM has been initialized. Be sure to call registration.kill_jvm() or slide_io.kill_jvm() at the end of your script."?

Do I need to do something to fix them?

Thanks

Marco

[marco1108]

Hi Chandler,

Another (hopefully last) question: I seem to be unable to modify the new .geojson once imported on QuPath. As if the annotations were blocked - but they aren't. Do you know why this could be?

Many thanks

M

[cdgatenbee]

Hi @marco1108,
Sorry, I still haven't been able to work on finding the fix about the SLF4J messages, but I don't think it should affect the performance of the registration. The last line is just a separate message from valis to provide a reminder about the need to shut down the jvm, because if it's not shut down the script won't finish. So, nothing really to do about those messages, everything should still work ok. About your other question, I'm sorry to say that I haven't done much in QuPath, but using the example registered geojson I have, I wasn't able to reproduce your error. The imported annotations appeared unlocked, and I was able to move the corners of simple polygons (e.g. squares, triangles, etc...), and drag around some of the ones that looked hand drawn (QuPath version 0.4.2). I should mention that this was using the registered annotations produced by valis version 1.0.0rc15, which has the Slide.warp_geojson_from_to method, which might be a bit different than the code you are using above (I realized this method doesn't show up on readthedocs, so I'll have to fix that). Not sure if that would resolve this issue, though. Sorry I couldn't be of more help here, but please do let me know if you are able to figure out how to modify your annotations.

Best,
-Chandler

[marco1108]

Hi Chandler,
Found out what was happening, again nothing to do with the script. I found out that if the annotations are locked in the original file from which you extract the .geojson, then you are unable to modify them on the new .geojson, even though they appear unlocked. Instead if I import them as unlock annotations everything works fine.
Thanks for helping! And thanks for creating this amazing script, it will save endless hours of work to our lab!
M