readme.md

Purpose

Export ROIs generated in the PathViewer UI for downstream analysis.

Execution

1. In google chrome, go the hms omero and make sure you are logged in
2. In any of your chrome tab/window, hit F12 to launch the developer tools
3. In the console tab of the DevTools paste in the contents of the export_omero_roi.js script
4. Replace the 7-digit omero image ID in the imgIds list with your image ID
5. Hit enter to run, it will ask you to save a csv file, your ROIs will be in there.

Output

A csv file image_name-image_id-rois.csv, each row is one ROI and the column headers are -

Id	Name	Text	type	all_points	X	Y	RadiusX	RadiusY	Width	Height	all_transforms

• Id: id of the ROI
• Name: name of the ROI
• Text: the text label from the user when an ROI is created
• type: type of the ROI, refer to reference 1
• all_points: coordinates for all corners "x1,y1 x2,y2 x3,y3 ..."
• X, Y, RadiusX, RadiusY, Width, Height: refer to reference 1
• all_transforms: refer to reference 2, "A00,A01,A02,A10,A11,A12,0,0,1"
• If the property does not exist for the shape, -1 is used as placeholder
• For ellipse, the four points are the vertices and co-vertices

Reference

1. OMERO ROI model
2. Affine Transformations of ROI Shapes

export_omero_roi.js

var imgIds = [
    1203716, 1203756, 1203720,
];
const omeroUrls = {
    hms: 'https://omero.hms.harvard.edu',
    idp: 'https://idp.tissue-atlas.org'
};
// choose omero instance between the HMS and the IDP OMERO
const whichOmero = 'hms';

// add 1 sec pause for each 10-download batch
imgIds.forEach((el, idx) => setTimeout(() => export_roi_by_id(el, omeroUrls[whichOmero.toLowerCase()]), idx * 200));

function export_roi_by_id(imgId, omeroUrl) {

    var omeroUrl = omeroUrl;
    var headers = ['Id', 'Name', 'Text', 'type', 'all_points', 'X', 'Y', 'RadiusX', 'RadiusY', 'Width', 'Height', 'all_transforms'];
    var url = `${omeroUrl}/api/v0/m/images/${imgId}/`;
    var urlRoi = `${omeroUrl}/api/v0/m/images/${imgId}/rois/?limit=500`;
    var imgName;

    fetch(url)
        .then(res => res.json())
        .then(resJson =>
            imgName = resJson.data.Name
        )
        .then(() => fetch(urlRoi)
            .then(res => res.json())
            .then(resJson =>
                resJson.data
                    .map(data => {
                        data.shapes[0].Name = data.Name || 'undefined';
                        data.shapes[0].Id = data['@id'];
                        return data.shapes;
                    })
                    .reduce((acc, val) => acc.concat(val), [])
                    // .filter(shape => shape['@type'].includes('Rect'))
                    .map(shape => {
                        shape.type = shape['@type'].split('#').pop();
                        shape.all_points = getPointsOfShape(shape);
                        shape.all_transforms = shape.Transform
                            ? [
                                shape.Transform.A00, shape.Transform.A01, shape.Transform.A02,
                                shape.Transform.A10, shape.Transform.A11, shape.Transform.A12,
                                0, 0, 1
                            ].join(',')
                            : -1;
                        return shape;
                    })
                    .map(shape => headers.map(header => JSON.stringify(shape[header]) || -1))
            )
            .then(rois => rois.map(e => e.join(",")).join("\n"))
            .then(roisStr =>
                "data:text/csv;charset=utf-8," +
                headers.join(',') + '\n' +
                roisStr.replace(/#/g, '-'))
            .then(fullStr => encodeURI(fullStr))
            .then(uri => {
                let link = document.createElement("a");
                link.setAttribute("href", uri);
                link.setAttribute("download", `${imgName}-${imgId}-rois.csv`);
                document.body.appendChild(link);
                link.click();
            })
        )
}

function getPointsOfShape(shape) {
    let all_points = [];
    if (shape['@type'].includes('Poly')) {
        const T = shape.Transform;
        all_points = shape.Points.split(' ').map(
            xy => xy.split(',').map(cc => parseFloat(cc))
        ).map(xy => transformPoint(...xy, T));
    }
    if (shape['@type'].includes('Rect')) {
        const [X, Y] = [shape.X, shape.Y];
        const [W, H] = [shape.Width, shape.Height];
        const T = shape.Transform;
        all_points = [
            transformPoint(X, Y, T),
            transformPoint(X + W, Y, T),
            transformPoint(X + W, Y + H, T),
            transformPoint(X, Y + H, T)
        ];
    }
    if (shape['@type'].includes('Ellipse')) {
        const [Rx, Ry] = [shape.RadiusX, shape.RadiusY];
        const [Cx, Cy] = [shape.X, shape.Y];
        const T = shape.Transform;
        all_points = [
            transformPoint(Cx + Rx, Cy, T),
            transformPoint(Cx - Rx, Cy, T),
            transformPoint(Cx, Cy + Ry, T),
            transformPoint(Cx, Cy - Ry, T),
        ];
    }
    if (shape['@type'].includes('Point')) {
        const T = shape.Transform;
        all_points = [
            transformPoint(shape.X, shape.Y, T)
        ];
    }
    if (!all_points.length) {
        console.warn(
            `ROI ${shape.Name} of type ${shape.type} is not fully supported yet in this script, validation is needed.`
        );
    }
    all_points = all_points.map(point => point.join(',')).join(' ');
    return all_points;
}
function transformPoint(x, y, t) {
    t = t
        ? t
        : { A00: 1, A01: 0, A02: 0, A10: 0, A11: 1, A12: 0 };
    const svg = document.createElementNS("http://www.w3.org/2000/svg", "svg");
    const point = svg.createSVGPoint();
    point.x = x;
    point.y = y;
    const transform = svg.createSVGMatrix();
    transform.a = t.A00;
    transform.b = t.A10;
    transform.c = t.A01;
    transform.d = t.A11;
    transform.e = t.A02;
    transform.f = t.A12;
    return [
        point.matrixTransform(transform).x.toFixed(3),
        point.matrixTransform(transform).y.toFixed(3)
    ];
}

function getUpperLeft() {
    let args = [...arguments];
    return [
        Math.min(...args.map(el => el.x)),
        Math.min(...args.map(el => el.y))
    ];
}
function getLowerRight() {
    let args = [...arguments];
    return [
        Math.max(...args.map(el => el.x)),
        Math.max(...args.map(el => el.y))
    ];
}

omero_roi.py

import pandas as pd
import re
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import scipy.spatial.distance as sdistance


def parse_roi_points(all_points):
    return np.array(
        re.findall(r'-?\d+\.?\d+', all_points), dtype=float
    ).reshape(-1, 2)

def ellipse_points_to_patch(
    vertex_1, vertex_2,
    co_vertex_1, co_vertex_2,
    patch_kwargs={}
):
    """
    Parameters
    ----------
    vertex_1, vertex_2, co_vertex_1, co_vertex_2: array like, in the form of (x-coordinate, y-coordinate)

    """
    v_and_co_v = np.array([
        vertex_1, vertex_2,
        co_vertex_1, co_vertex_2
    ])
    centers = v_and_co_v.mean(axis=0)

    d = sdistance.cdist(v_and_co_v, v_and_co_v, metric='euclidean')
    width = d[0, 1]
    height = d[2, 3]

    vector_2 = v_and_co_v[1] - v_and_co_v[0]
    vector_2 /= np.linalg.norm(vector_2)

    angle = np.degrees(np.arccos([1, 0] @ vector_2))

    ellipse_patch = mpatches.Ellipse(
        centers, width=width, height=height, angle=angle,
        **patch_kwargs
    )
    return ellipse_patch

def add_mpatch(roi, patch_kwargs={}):
    roi = roi.copy()
    points = parse_roi_points(roi['all_points'])
    roi.loc['parsed_points'] = points

    roi_type = roi['type']
    if roi_type in ['Point', 'Line']:
        roi_mpatch = mpatches.Polygon(points, closed=False, **patch_kwargs)
    elif roi_type in ['Rectangle', 'Polygon', 'Polyline']:
        roi_mpatch = mpatches.Polygon(points, closed=True, **patch_kwargs)
    elif roi_type == 'Ellipse':
        roi_mpatch = ellipse_points_to_patch(*points, patch_kwargs=patch_kwargs)
    else:
        raise ValueError

    roi.loc['mpatch'] = roi_mpatch
    return roi

# use dask for parallel computing when num of points is large
import dask.dataframe as dd
import dask.diagnostics
def compute_inside_mask(roi, df, parallel=False):
    sub_df = df[['X_centroid', 'Y_centroid']][
        df.X_centroid.between(roi.x_min, roi.x_max, inclusive='both') * 
        df.Y_centroid.between(roi.y_min, roi.y_max, inclusive='both')
    ]
    if parallel & (sub_df.shape[0] > 100000):
        dd_sub_df = dd.from_pandas(sub_df, chunksize=10000)
        with dask.diagnostics.ProgressBar():
            inside = dd_sub_df.map_partitions(
                roi.mpatch.contains_points, meta=(None, bool)
            ).compute(scheduler='processes')
    else:
        inside = roi.mpatch.contains_points(sub_df)
    mask = df.X_centroid > df.X_centroid + 1
    mask.loc[sub_df[inside].index] = True
    return mask.values

def demo():
    # load single-cell quantifacation table
    ss = pd.read_csv('quantification/cellRingMask/Lung_061_CFX16a_Full_C11@20201024_195941_804220.csv', index_col=0)

    # load ROI table exported from OMERO
    _roi_table = pd.read_csv('Lung_061_CFX16a_Full_C11@20201024_195941_804220.ome.tiff-1327440-rois.csv')

    # set patch properties for visualization
    patch_kwargs = dict(edgecolor='red', facecolor='none')

    # add 'mpatch' column to the returned table
    roi_table = _roi_table.apply(add_mpatch, axis=1, patch_kwargs=patch_kwargs)

    # precompute bounding box for better performance
    roi_table.loc[:, ['x_min', 'y_min', 'x_max', 'y_max']] = np.hstack([
        [*roi_table['parsed_points'].apply(np.min, axis=0)],
        [*roi_table['parsed_points'].apply(np.max, axis=0)]
    ])

    # add 'inside_mask' column to the table
    roi_table.loc[:, 'inside_mask'] = roi_table.apply(
        compute_inside_mask, df=ss, parallel=True, axis=1
    )
    inside_mask_all = np.logical_or.reduce(roi_table['inside_mask'].to_list())

    plt.figure()
    plt.scatter(
        ss['X_centroid'], ss['Y_centroid'], 
        marker='.', linewidths=0, s=2, c=inside_mask_all
    )
    plt.gca().set_aspect('equal')
    plt.gca().invert_yaxis()
    for p in roi_table['mpatch']:
        plt.gca().add_patch(p)

omero_roi_to_idp.py

from omero.gateway import BlitzGateway
import omero
from omero.model.enums import UnitsLength


def clean_shape(shape):
    default_vals = {
        '_details': omero.model.DetailsI(),
        '_id': None,
        '_roi': None,
    }
    def reset(target):
        for k, v in default_vals.items():
            if hasattr(target, k):
                setattr(target, k, v)
    reset(shape)
    for key in shape.__dict__.keys():
        secondary_entry = getattr(shape, key)
        reset(secondary_entry)
    return shape


def change_session_group_to_img_id(conn, img_id):
    conn.SERVICE_OPTS.setOmeroGroup('-1')

    image = conn.getObject("Image", int(img_id))
    group_id = image.getDetails().getGroup().getId()
    # Switch session group to image's group 
    conn.setGroupForSession(group_id)
    return 


HMS_SESSION_ID = '72d725ae-a082-4101-b6bf-1e87a4470383'
# Create a connection
# ===================
conn_hms = BlitzGateway(
    HMS_SESSION_ID, HMS_SESSION_ID,
    host='omero-app.hms.harvard.edu',
    port=4064
)
conn_hms.connect()


IDP_SESSION_ID = '4cd42002-9915-401d-9ad5-7c3eff21b57b'
# Create a connection
# ===================
conn_idp = BlitzGateway(
    IDP_SESSION_ID, IDP_SESSION_ID,
    host='idp.tissue-atlas.org',
)
conn_idp.connect()


HMS_IMAGE_ID = 1216351
IDP_IMAGE_ID = 53

# query source rois
change_session_group_to_img_id(conn_hms, HMS_IMAGE_ID)
roi_service = conn_hms.getRoiService()
result = roi_service.findByImage(HMS_IMAGE_ID, None)


roi_names = [r.name for r in result.rois]
roi_shapes = [r.copyShapes() for r in result.rois]
print(f'# of rois from source ({HMS_IMAGE_ID}): {len(result.rois)}')
      
# discard elements related to user/session info from source rois
roi_shapes = [
    [clean_shape(s) for s in shapes]
    for shapes in roi_shapes
]

# create new rois for target image from source rois
rois_for_target = []
for n, ss in zip(roi_names, roi_shapes):
    roi_new = omero.model.RoiI()
    roi_new.name = n
    for s in ss:
        roi_new.addShape(s)
    rois_for_target.append(roi_new)

print(f'# of rois to target ({IDP_IMAGE_ID}): {len(rois_for_target)}')


# update service @ target host and helper func to add rois
updateService = conn_idp.getUpdateService()
def create_roi(img, rois):
    for roi in rois:
        roi.setImage(img._obj)
        updateService.saveAndReturnObject(roi)
    return


# add the rois to target image
change_session_group_to_img_id(conn_idp, IDP_IMAGE_ID)
image = conn_idp.getObject("Image", int(IDP_IMAGE_ID))
create_roi(image, rois_for_target)

conn_hms.close()
conn_idp.close()

roi2mask.py

from matplotlib.pyplot import fill
import omero_roi
import pandas as pd
import numpy as np
import cv2


def roi_path_to_mask(
    roi_path,
    mask_shape=None,
    fill_value=None
):
    _roi = pd.read_csv(roi_path)
    roi = _roi.apply(omero_roi.add_mpatch, axis=1)

    if mask_shape is None:
        mask_shape = [
            roi.parsed_points.apply(
                lambda x: np.array(x).max(axis=0)[i]
            ).max()
            for i in range(2)
        ][::-1]
    h, w = mask_shape
    
    if fill_value is None:
        fill_value = 1
    fill_value = int(fill_value)

    mask = np.zeros((int(h), int(w)), np.uint8)
    cv2.fillPoly(
        mask,
        roi.parsed_points.apply(lambda x: np.round(x).astype(int)),
        fill_value
    )
    return mask

validation.py

import pandas as pd
import re
import skimage.io
from skimage.draw import polygon
import numpy as np


img = skimage.io.imread('16_c0.tif')
table = pd.read_csv('16.ome.tif-1207098-rois.csv')

all_points = [
    np.array(re.findall(r'\d+\.\d+', s))
        .astype(np.float64)
        .reshape(-1, 2)
    for s in table['all_points']
]

def find_bound(coors):
    # return [[x_min, ymin], [x_max, y_max]]
    return (
        coors.min(axis=0).astype(np.int64),
        coors.max(axis=0).astype(np.int64)
        )

bounds = [find_bound(i) for i in all_points]

def crop_roi(idx):
    points = all_points[idx]
    bound = bounds[idx]

    # Add 1 to upper bounds for subpixels
    roi = img[
        bound[0][1]:bound[1][1] + 1,
        bound[0][0]:bound[1][0] + 1
    ]

    mask = np.zeros_like(roi)
    rr, cc = polygon((points - bound[0])[..., 1], (points - bound[0])[..., 0])
    mask[rr, cc] = 1
    
    skimage.io.imsave(
        '{:03}-{}.tif'.format(idx + 1, table['Name'][idx]),
        roi * mask
    )

for i in range(len(table['Name'])):
    crop_roi(i)