diff --git a/transforms/images/polus-image-registration-plugin/.bumpversion.cfg b/transforms/images/polus-image-registration-plugin/.bumpversion.cfg new file mode 100644 index 000000000..1dd0a7024 --- /dev/null +++ b/transforms/images/polus-image-registration-plugin/.bumpversion.cfg @@ -0,0 +1,16 @@ +[bumpversion] +current_version = 0.3.6 +commit = False +tag = False + +[bumpversion:file:VERSION] + +[bumpversion:file:pyproject.toml] +search = version = "{current_version}" +replace = version = "{new_version}" + +[bumpversion:file:plugin.json] + +[bumpversion:file:ict.yaml] + +[bumpversion:file:imageregistration.cwl] diff --git a/transforms/images/polus-image-registration-plugin/Dockerfile b/transforms/images/polus-image-registration-plugin/Dockerfile index 3f94d8778..09aca4f24 100644 --- a/transforms/images/polus-image-registration-plugin/Dockerfile +++ b/transforms/images/polus-image-registration-plugin/Dockerfile @@ -1,21 +1,16 @@ - -FROM polusai/bfio:2.1.9 - -ARG EXEC_DIR="/opt/executables" -ARG DATA_DIR="/data" - -RUN mkdir -p ${EXEC_DIR} \ - && mkdir -p ${DATA_DIR}/inputs \ - && mkdir ${DATA_DIR}/outputs - -RUN apt-get update && \ - apt-get -y install libglib2.0-0 - -COPY src ${EXEC_DIR}/ +# Build from repo root (monorepo) or from this tool directory — both work. +FROM polusai/bfio:2.5.0 +ENV EXEC_DIR="/opt/executables" POLUS_IMG_EXT=".ome.tif" POLUS_TAB_EXT=".csv" POLUS_LOG="INFO" WORKDIR ${EXEC_DIR} - -RUN pip3 install -r ${EXEC_DIR}/requirements.txt --no-cache-dir - -COPY VERSION / - -ENTRYPOINT ["python3", "/opt/executables/main.py"] \ No newline at end of file +ENV TOOL_DIR="transforms/images/polus-image-registration-plugin" +RUN mkdir -p image-tools +COPY . ${EXEC_DIR}/image-tools +RUN apt-get update && apt-get install -y --no-install-recommends libglib2.0-0 \ + && rm -rf /var/lib/apt/lists/* +RUN pip3 install -U pip setuptools wheel \ + && python3 -c 'import sys; assert sys.version_info>=(3,11)' \ + && R="${EXEC_DIR}/image-tools" && M="$R/$TOOL_DIR" \ + && if [ -f "$M/pyproject.toml" ]; then pip3 install --no-cache-dir "$M"; \ + else pip3 install --no-cache-dir "$R"; fi +ENTRYPOINT ["python3", "-m", "main"] +CMD ["--help"] diff --git a/transforms/images/polus-image-registration-plugin/README.md b/transforms/images/polus-image-registration-plugin/README.md index 4aac366aa..14fa6710e 100644 --- a/transforms/images/polus-image-registration-plugin/README.md +++ b/transforms/images/polus-image-registration-plugin/README.md @@ -1,6 +1,6 @@ # Polus Projective Transformation Image Registration Plugin -WIPP Plugin Title : Image Registration Plugin +WIPP Plugin Title : Image Registration Plugin Contact [Gauhar Bains](mailto:gauhar.bains@labshare.org) for more information. @@ -8,53 +8,53 @@ For more information on WIPP, visit the [official WIPP page](https://isg.nist.go ## Description -This plugin registers an image collection. First it parses the image collection using `parser.py` into registration sets. Each registration set consists of: moving image, template image, similiar transformation images. The registration algorithm(explained in the next section, can be found in `image_registration.py`) registers the moving image with template image and stores the transformation required to do so. This stored transformation is used to transform each image in similar transformation list. +This plugin registers an image collection. First it parses the image collection using `parser.py` into registration sets. Each registration set consists of: moving image, template image, similiar transformation images. The registration algorithm(explained in the next section, can be found in `image_registration.py`) registers the moving image with template image and stores the transformation required to do so. This stored transformation is used to transform each image in similar transformation list. ## Algorithm ### To do -1. Find a better way to handle poorly correlated transforms. +1. Find a better way to handle poorly correlated transforms. 2. Find a more scalable approach to do rough transformation. The algorithm has been tested on images of size around 1.2 Gigapixel but a better approach may be needed for images significantly larger than these. -### Parsing -The parsing algorithm uses the functions from the `file_pattern utility`. It takes the following inputs : Filename pattern, registration variable, similar transformation variable. The registration variable helps determine the moving and the template images where as the similar transformation variable helps determine the similar transformation images. Note: The code produces the expected output when len(registration_variable)==len(similarity_variable)==1. The code will NOT spit out an error when the more than one variable is passed as registration or similarity variable, but additional testing needs to be done for this usecase. +### Parsing +The parsing algorithm uses the functions from the `file_pattern utility`. It takes the following inputs : Filename pattern, registration variable, similar transformation variable. The registration variable helps determine the moving and the template images where as the similar transformation variable helps determine the similar transformation images. Note: The code produces the expected output when len(registration_variable)==len(similarity_variable)==1. The code will NOT spit out an error when the more than one variable is passed as registration or similarity variable, but additional testing needs to be done for this usecase. -Some sample text files can be found in the examples folder. Short example shown below: +Some sample text files can be found in the examples folder. Short example shown below: -Parsing example : - -`Inputs:` -Filepattern : `x{xxx}_y{yyy}_z{zzz}_c{ccc}_t{ttt}.ome.tif` -Registration_variable : `t` -similar_transformation_variable : `c` -template_ : `x001_y001_z001_c001_t001.ome.tif` +Parsing example : -`Output set 1 :` -Template Image: x001_y001_z001_c001_t001.ome.tif -Moving Image: x001_y001_z001_c001_t002.ome.tif -Similar Transformation Images : [ x001_y001_z001_c002_t002.ome.tif , x001_y001_z001_c003_t002.ome.tif ] +`Inputs:` +Filepattern : `x{xxx}_y{yyy}_z{zzz}_c{ccc}_t{ttt}.ome.tif` +Registration_variable : `t` +similar_transformation_variable : `c` +template_ : `x001_y001_z001_c001_t001.ome.tif` -`Output set 2:` -Template Image: x001_y002_z001_c001_t001.ome.tif -Moving Image: x001_y002_z001_c001_t002.ome.tif -Similar Transformation Images : [ x001_y002_z001_c002_t002.ome.tif , x001_y002_z001_c003_t002.ome.tif ] +`Output set 1 :` +Template Image: x001_y001_z001_c001_t001.ome.tif +Moving Image: x001_y001_z001_c001_t002.ome.tif +Similar Transformation Images : [ x001_y001_z001_c002_t002.ome.tif , x001_y001_z001_c003_t002.ome.tif ] +`Output set 2:` +Template Image: x001_y002_z001_c001_t001.ome.tif +Moving Image: x001_y002_z001_c001_t002.ome.tif +Similar Transformation Images : [ x001_y002_z001_c002_t002.ome.tif , x001_y002_z001_c003_t002.ome.tif ] -### Registration -The registration algorithm is present in `image_registration.py`. It uses projective transformation(Homography matrix) to alter the moving image and align it with the reference image. Background information about homography can be found here: https://en.wikipedia.org/wiki/Homography . -The moving image undergoes 2 transformations: + +### Registration +The registration algorithm is present in `image_registration.py`. It uses projective transformation(Homography matrix) to alter the moving image and align it with the reference image. Background information about homography can be found here: https://en.wikipedia.org/wiki/Homography . +The moving image undergoes 2 transformations: 1. `Rough Transformation` : In this the whole moving image is transformed using the homography matrix calculated between the entire moving and template image. -2. `Fine Transformation` : To carry out fine transformation ,the homography matrix is found between the corresponding tiles of the roughly transformed moving image and the template image. Each image is divided into 4 tiles. +2. `Fine Transformation` : To carry out fine transformation ,the homography matrix is found between the corresponding tiles of the roughly transformed moving image and the template image. Each image is divided into 4 tiles. + +To find the homography matrix(for fine or rough tranformation), we need coordinates of atleast 4 matching points in the template and the moving image. To do this the ORB feature detector has been used. However, its computationally very expensive to run feature matching on large images(our test data consists of 1.3 gigapixel images). To overcome this, the homography matrix at every step of our algorithm has been calculated between scaled down versions( 16 * 16 times smaller) of the respective images. To use this homography matrix on actual sized images, the matrix is scaled up using a scale matrix. The proof for upscaling a homography matrix is shown below. -To find the homography matrix(for fine or rough tranformation), we need coordinates of atleast 4 matching points in the template and the moving image. To do this the ORB feature detector has been used. However, its computationally very expensive to run feature matching on large images(our test data consists of 1.3 gigapixel images). To overcome this, the homography matrix at every step of our algorithm has been calculated between scaled down versions( 16 * 16 times smaller) of the respective images. To use this homography matrix on actual sized images, the matrix is scaled up using a scale matrix. The proof for upscaling a homography matrix is shown below. +`Proof` : -`Proof` : +Credit for proof : https://stackoverflow.com/questions/21019338/how-to-change-the-homography-with-the-scale-of-the-image/56623249 + +![homography](https://user-images.githubusercontent.com/48079888/78402511-b04d8200-75c8-11ea-9d22-cee13f3912db.gif) -Credit for proof : https://stackoverflow.com/questions/21019338/how-to-change-the-homography-with-the-scale-of-the-image/56623249 - -![homography](https://user-images.githubusercontent.com/48079888/78402511-b04d8200-75c8-11ea-9d22-cee13f3912db.gif) - @@ -83,4 +83,3 @@ This plugin takes one input argument and one output argument: | `--template` | Template image to be used for image registration | Input | string | | `--TransformationVariable` | variable to help identify which images have similar transformation | Input | string | | `--outDir` | Output collection | Output | collection | - diff --git a/transforms/images/polus-image-registration-plugin/VERSION b/transforms/images/polus-image-registration-plugin/VERSION index c2c0004f0..449d7e73a 100644 --- a/transforms/images/polus-image-registration-plugin/VERSION +++ b/transforms/images/polus-image-registration-plugin/VERSION @@ -1 +1 @@ -0.3.5 +0.3.6 diff --git a/transforms/images/polus-image-registration-plugin/build-docker.sh b/transforms/images/polus-image-registration-plugin/build-docker.sh index 0809b7e34..43dd281d0 100755 --- a/transforms/images/polus-image-registration-plugin/build-docker.sh +++ b/transforms/images/polus-image-registration-plugin/build-docker.sh @@ -1,4 +1,4 @@ #!/bin/bash version=$(=1.26.4", + "opencv-python-headless>=4.7.0", + "filepattern>=2.2.1", +] + +[build-system] +requires = ["hatchling"] +build-backend = "hatchling.build" + +[tool.hatch.build.targets.wheel] +packages = ["src"] diff --git a/transforms/images/polus-image-registration-plugin/src/__init__.py b/transforms/images/polus-image-registration-plugin/src/__init__.py new file mode 100644 index 000000000..c197bd65d --- /dev/null +++ b/transforms/images/polus-image-registration-plugin/src/__init__.py @@ -0,0 +1 @@ +"""Image registration plugin.""" diff --git a/transforms/images/polus-image-registration-plugin/src/image_registration.py b/transforms/images/polus-image-registration-plugin/src/image_registration.py index fdc353ba2..b775f331c 100644 --- a/transforms/images/polus-image-registration-plugin/src/image_registration.py +++ b/transforms/images/polus-image-registration-plugin/src/image_registration.py @@ -1,440 +1,670 @@ -import cv2, argparse, logging -import numpy as np -from bfio.bfio import BioReader, BioWriter -from pathlib import Path +"""Image registration: ORB/FLANN matching, homography, and tiled bfio read/write.""" +from __future__ import annotations + +import argparse +import logging from concurrent.futures import ThreadPoolExecutor from multiprocessing import cpu_count +from pathlib import Path +from typing import cast + +import cv2 +import numpy as np +from bfio.bfio import BioReader +from bfio.bfio import BioWriter -def corr2(a,b): - """corr2 Calculate correlation between 2 images +CORRELATION_FALLBACK_THRESHOLD = 0.4 +TILE_MERGE_THRESHOLD = 3072 +TARGET_PIXELS = 5_000_000 # 5 megapixels downscale target - Inputs: - a (np.ndarray): A 2-dimensional numpy array - b (np.ndarray): A 2-dimensional numpy array - Outputs: - float: the correlation between a and b +def corr2(a: np.ndarray, b: np.ndarray) -> float: + """Correlation between two same-shaped 2D arrays (Pearson on flattened data). + + Args: + a: First image or tile. + b: Second image or tile. + + Returns: + Scalar correlation coefficient. """ - - c = np.sum(a)/np.size(a) - d = np.sum(b)/np.size(b) - + c = np.sum(a) / np.size(a) + d = np.sum(b) / np.size(b) + c = a - c d = b - d - - r = (c*d).sum() / np.sqrt((c*c).sum() * (d*d).sum()) - - return r - -def get_transform(moving_image,reference_image,max_val,min_val,method): - """get_transform Calculate homography matrix transform - - This function registers the moving image with reference image - - Inputs: - moving_image = Image to be transformed - reference_image= reference Image - Outputs: - homography= transformation applied to the moving image + + return float((c * d).sum() / np.sqrt((c * c).sum() * (d * d).sum())) + + +def get_transform( + moving_image: np.ndarray, + reference_image: np.ndarray, + max_val: float, + min_val: float, + method: str, +) -> np.ndarray | None: + """Estimate homography or affine map from moving_image to reference_image. + + Args: + moving_image: Image to align. + reference_image: Fixed reference. + max_val: Normalization upper bound. + min_val: Normalization lower bound. + method: One of ``Projective``, ``Affine``, or ``PartialAffine``. + + Returns: + Transform matrix, or ``None`` if matching fails. """ - # height, width of the reference image - height, width = reference_image.shape # max number of features to be calculated using ORB - max_features=500000 + max_features = 500000 # initialize orb feature matcher orb = cv2.ORB_create(max_features) - + # Normalize images and convert to appropriate type - moving_image_norm = cv2.GaussianBlur(moving_image,(3,3),0) - reference_image_norm = cv2.GaussianBlur(reference_image,(3,3),0) - moving_image_norm = (moving_image_norm-min_val)/(max_val-min_val) + moving_image_norm = cv2.GaussianBlur(moving_image, (3, 3), 0) + reference_image_norm = cv2.GaussianBlur(reference_image, (3, 3), 0) + moving_image_norm = (moving_image_norm - min_val) / (max_val - min_val) moving_image_norm = (moving_image_norm * 255).astype(np.uint8) - reference_image_norm = (reference_image_norm-min_val)/(max_val-min_val) + reference_image_norm = (reference_image_norm - min_val) / (max_val - min_val) reference_image_norm = (reference_image_norm * 255).astype(np.uint8) - + # find keypoints and descriptors in moving and reference image keypoints1, descriptors1 = orb.detectAndCompute(moving_image_norm, None) keypoints2, descriptors2 = orb.detectAndCompute(reference_image_norm, None) - + # Escape if one image does not have descriptors - if not (isinstance(descriptors1,np.ndarray) and isinstance(descriptors2,np.ndarray)): + d1_ok = isinstance(descriptors1, np.ndarray) + d2_ok = isinstance(descriptors2, np.ndarray) + desc_ok = d1_ok and d2_ok + if not desc_ok: return None - + # match and sort the descriptos using hamming distance - flann_params= dict(algorithm = 6, # FLANN_INDEX_LSH - table_number = 6, - key_size = 12, - multi_probe_level = 1) + flann_params = { + "algorithm": 6, # FLANN_INDEX_LSH + "table_number": 6, + "key_size": 12, + "multi_probe_level": 1, + } matcher = cv2.FlannBasedMatcher(flann_params, {}) matches = matcher.match(descriptors1, descriptors2, None) matches.sort(key=lambda x: x.distance, reverse=False) - + # extract top 25% of matches - good_match_percent=0.25 - numGoodMatches = int(len(matches) * good_match_percent) - matches = matches[:numGoodMatches] - + good_match_percent = 0.25 + num_good_matches = int(len(matches) * good_match_percent) + matches = matches[:num_good_matches] + # extract the point coordinates from the keypoints points1 = np.zeros((len(matches), 2), dtype=np.float32) - points2 = np.zeros((len(matches), 2), dtype=np.float32) + points2 = np.zeros((len(matches), 2), dtype=np.float32) for i, match in enumerate(matches): points1[i, :] = keypoints1[match.queryIdx].pt - points2[i, :] = keypoints2[match.trainIdx].pt - + points2[i, :] = keypoints2[match.trainIdx].pt + # If no matching points, return None - if points1.shape[0]==0 or points2.shape[0]==0: + if points1.shape[0] == 0 or points2.shape[0] == 0: return None - + # calculate the homography matrix - if method=='Projective': + homography: np.ndarray + if method == "Projective": homography, _ = cv2.findHomography(points1, points2, cv2.RANSAC) - elif method=='Affine': + elif method == "Affine": homography, _ = cv2.estimateAffine2D(points1, points2, method=cv2.RANSAC) - elif method=='PartialAffine': + elif method == "PartialAffine": homography, _ = cv2.estimateAffinePartial2D(points1, points2, method=cv2.RANSAC) - + else: + return None + return homography -def get_scale_factor(height,width): - """ - This function returns the appropriate scale factor w.r.t to - a target size. Target size has been fixed to 5 megapixels. - - Inputs: - height (int): Image height - width (int): Image width - Outputs: - scale factor - """ - TARGET_SIZE=5000000 # 5 megapixels - scale_factor=np.sqrt((height*width)/TARGET_SIZE) - return int(scale_factor) if scale_factor>1 else 1 - -def get_scaled_down_images(image,scale_factor,get_max=False): + +def get_scale_factor(height: int, width: int) -> int: + """Scale factor to cap effective resolution near ``TARGET_PIXELS``. + + Args: + height: Image height. + width: Image width. + + Returns: + Integer scale factor (at least 1). """ - This function returns the scaled down version of an image. - Inputs: - image : A BioReader object - scale_factor : the factor by which the image needs - to be scaled down - Outputs: - rescaled_image: scaled down version of the input image + target_size = TARGET_PIXELS + scale_factor = np.sqrt((height * width) / target_size) + return int(scale_factor) if scale_factor > 1 else 1 + + +def get_scaled_down_images( + image: BioReader, + scale_factor: int, + get_max: bool = False, +) -> np.ndarray | tuple[np.ndarray, float, float]: + """Return a downscaled view of a large BioReader image (optionally min/max). + + Args: + image: Open ``BioReader`` instance. + scale_factor: Integer downscale factor. + get_max: If True, also return max and min pixel values from tiles. + + Returns: + Downscaled array, or ``(array, max_val, min_val)`` when ``get_max`` is True. """ - # Calculate scaling variables - stride = int(scale_factor * np.floor(4096/scale_factor)) - width = np.ceil(image.num_y()/scale_factor).astype(int) - height = np.ceil(image.num_x()/scale_factor).astype(int) - + stride = int(scale_factor * np.floor(4096 / scale_factor)) + width = np.ceil(image.num_y() / scale_factor).astype(int) + height = np.ceil(image.num_x() / scale_factor).astype(int) + # Initialize the output - rescaled_image = np.zeros((width,height),dtype=image._pix['type']) - - # If max value is requested, initialize the variables - if get_max: - max_val = 0 - min_val = np.inf - - def load_and_scale(X,Y,x,y,get_max=get_max,reader=image,scale_factor=scale_factor,rescaled_image=rescaled_image): - """load_and_scale Load a section of an image and downscale - - This is a transient method, and only works within the get - scaled_down_images method. - It's used to thread out loading and downscaling of large images. - - """ - + rescaled_image = np.zeros((width, height), dtype=image._pix["type"]) + + def load_and_scale( # noqa: PLR0913 + x_read_bounds: list[int], + y_read_bounds: list[int], + x_write_bounds: list[int], + y_write_bounds: list[int], + get_max: bool = get_max, + reader: BioReader = image, + scale_factor: int = scale_factor, + rescaled_image: np.ndarray = rescaled_image, + ) -> tuple[float, float] | None: + """Load a tile, downscale, and write into ``rescaled_image`` (thread worker).""" # Read an image tile - tile = reader.read_image(X=X,Y=Y,Z=[0,1],C=[0],T=[0]).squeeze() - + tile = reader.read_image( + X=x_read_bounds, + Y=y_read_bounds, + Z=[0, 1], + C=[0], + T=[0], + ).squeeze() + # Average the image for scaling - blurred_image = cv2.boxFilter(tile,-1,(scale_factor,scale_factor)) - + blurred_image = cv2.boxFilter(tile, -1, (scale_factor, scale_factor)) + # Collect pixels for downscaled image - rescaled_image[y[0]:y[1],x[0]:x[1]] = blurred_image[::scale_factor,::scale_factor] - + y0, y1 = y_write_bounds[0], y_write_bounds[1] + x0, x1 = x_write_bounds[0], x_write_bounds[1] + rescaled_image[y0:y1, x0:x1] = blurred_image[::scale_factor, ::scale_factor] + if get_max: - return np.max(tile),np.min(tile) - else: - return None - + return np.max(tile), np.min(tile) + return None + # Load and downscale the image threads = [] - with ThreadPoolExecutor(max([cpu_count()//2,1])) as executor: - for x in range(0,image.num_x(),stride): - x_max = np.min([x+stride,image.num_x()]) # max x to load - xi = int(x//scale_factor) # initial scaled x-index - xe = int(np.ceil(x_max/scale_factor)) # ending scaled x-index - for y in range(0,image.num_y(),stride): - y_max = np.min([y+stride,image.num_y()]) # max y to load - yi = int(y//scale_factor) # initial scaled y-index - ye = int(np.ceil(y_max/scale_factor)) # ending scaled y-index - - threads.append(executor.submit(load_and_scale,[x,x_max],[y,y_max],[xi,xe],[yi,ye])) - + with ThreadPoolExecutor(max([cpu_count() // 2, 1])) as executor: + for x in range(0, image.num_x(), stride): + x_max = np.min([x + stride, image.num_x()]) # max x to load + xi = int(x // scale_factor) # initial scaled x-index + xe = int(np.ceil(x_max / scale_factor)) # ending scaled x-index + for y in range(0, image.num_y(), stride): + y_max = np.min([y + stride, image.num_y()]) # max y to load + yi = int(y // scale_factor) # initial scaled y-index + ye = int(np.ceil(y_max / scale_factor)) # ending scaled y-index + + threads.append( + executor.submit( + load_and_scale, + [x, x_max], + [y, y_max], + [xi, xe], + [yi, ye], + ), + ) + # Return max and min values if requested if get_max: - results = [thread.result() for thread in threads] - max_val = max(result[0] for result in results) - min_val = min(result[1] for result in results) - return rescaled_image,max_val,min_val - else: - return rescaled_image + tile_ranges: list[tuple[float, float]] = [] + for thread in threads: + out = cast(tuple[float, float], thread.result()) + tile_ranges.append((float(out[0]), float(out[1]))) + max_val = max(t[0] for t in tile_ranges) + min_val = min(t[1] for t in tile_ranges) + return rescaled_image, max_val, min_val + return rescaled_image + -def register_image(br_ref,br_mov,bw,Xt,Yt,Xm,Ym,x,y,X_crop,Y_crop,max_val,min_val,method): - """register_image Register one section of two images +def register_image( # noqa: PLR0913 + br_ref: BioReader, + br_mov: BioReader, + bw: BioWriter, + xt_bounds: list[int], + yt_bounds: list[int], + xm_bounds: list[int], + ym_bounds: list[int], + tile_x: int, + tile_y: int, + x_crop_bounds: list[int], + y_crop_bounds: list[int], + max_val: float, + min_val: float, + method: str, + rough_homography_upscaled: np.ndarray, +) -> np.ndarray: + """Register one tile pair and write the warped moving tile to ``bw``. - This method is designed to be used within a thread. It registers - one section of two different images, saves the output, and - returns the homography matrix used to transform the image. + Args: + br_ref: Reference ``BioReader``. + br_mov: Moving ``BioReader``. + bw: Output ``BioWriter``. + xt_bounds: Reference X span ``[start, end)``. + yt_bounds: Reference Y span ``[start, end)``. + xm_bounds: Moving X span. + ym_bounds: Moving Y span. + tile_x: Output tile X index. + tile_y: Output tile Y index. + x_crop_bounds: Crop range in X after warp. + y_crop_bounds: Crop range in Y after warp. + max_val: Intensity normalization max. + min_val: Intensity normalization min. + method: Transform kind. + rough_homography_upscaled: Fallback transform if local match fails. + Returns: + Homography or affine matrix used for this tile. """ - # Load a section of the reference and moving images - ref_tile = br_ref.read_image(X=[Xt[0],Xt[1]],Y=[Yt[0],Yt[1]],Z=[0,1],C=[0],T=[0]).squeeze() - mov_tile = br_mov.read_image(X=[Xm[0],Xm[1]],Y=[Ym[0],Ym[1]],Z=[0,1],C=[0],T=[0]).squeeze() - + ref_tile = br_ref.read_image( + X=[xt_bounds[0], xt_bounds[1]], + Y=[yt_bounds[0], yt_bounds[1]], + Z=[0, 1], + C=[0], + T=[0], + ).squeeze() + mov_tile = br_mov.read_image( + X=[xm_bounds[0], xm_bounds[1]], + Y=[ym_bounds[0], ym_bounds[1]], + Z=[0, 1], + C=[0], + T=[0], + ).squeeze() + # Get the transformation matrix - projective_transform = get_transform(mov_tile,ref_tile,max_val,min_val,method) - + projective_transform = get_transform(mov_tile, ref_tile, max_val, min_val, method) + # Use the rough transformation matrix if no matrix was returned is_rough = False - if not isinstance(projective_transform,np.ndarray): + if not isinstance(projective_transform, np.ndarray): is_rough = True - projective_transform = Rough_Homography_Upscaled - + projective_transform = rough_homography_upscaled + # Transform the moving image - if method=='Projective': - transformed_image = cv2.warpPerspective(mov_tile,projective_transform,(Xt[1]-Xt[0],Yt[1]-Yt[0])) + if method == "Projective": + transformed_image = cv2.warpPerspective( + mov_tile, + projective_transform, + (xt_bounds[1] - xt_bounds[0], yt_bounds[1] - yt_bounds[0]), + ) else: - transformed_image = cv2.warpAffine(mov_tile,projective_transform,(Xt[1]-Xt[0],Yt[1]-Yt[0])) - + transformed_image = cv2.warpAffine( + mov_tile, + projective_transform, + (xt_bounds[1] - xt_bounds[0], yt_bounds[1] - yt_bounds[0]), + ) + # Determine the correlation between the reference and transformed moving image - corr = corr2(ref_tile,transformed_image) - + corr = corr2(ref_tile, transformed_image) + # If the correlation is bad, try using the rough transform instead - if corr < 0.4 and not is_rough: - if method=='Projective': - transformed_image = cv2.warpPerspective(mov_tile,Rough_Homography_Upscaled,(Xt[1]-Xt[0],Yt[1]-Yt[0])) + if corr < CORRELATION_FALLBACK_THRESHOLD and not is_rough: + if method == "Projective": + transformed_image = cv2.warpPerspective( + mov_tile, + rough_homography_upscaled, + (xt_bounds[1] - xt_bounds[0], yt_bounds[1] - yt_bounds[0]), + ) else: - transformed_image = cv2.warpAffine(mov_tile,Rough_Homography_Upscaled,(Xt[1]-Xt[0],Yt[1]-Yt[0])) - projective_transform = Rough_Homography_Upscaled - + transformed_image = cv2.warpAffine( + mov_tile, + rough_homography_upscaled, + (xt_bounds[1] - xt_bounds[0], yt_bounds[1] - yt_bounds[0]), + ) + projective_transform = rough_homography_upscaled + # Write the transformed moving image - bw.write_image(transformed_image[Y_crop[0]:Y_crop[1],X_crop[0]:X_crop[1],np.newaxis,np.newaxis,np.newaxis],X=[x],Y=[y]) - + bw.write_image( + transformed_image[ + y_crop_bounds[0] : y_crop_bounds[1], + x_crop_bounds[0] : x_crop_bounds[1], + np.newaxis, + np.newaxis, + np.newaxis, + ], + X=[tile_x], + Y=[tile_y], + ) + return projective_transform -def apply_transform(br_mov,bw,tiles,shape,transform,method): - """apply_transform Apply a transform to an image - This method is designed to be used within a thread. It loads - a section of an image, applies a transform, and saves the - transformed image to file. +def apply_transform( # noqa: PLR0913 + br_mov: BioReader, + bw: BioWriter, + tiles: tuple[list[int], list[int], list[int], list[int]], + shape: tuple[int, int, list[int], list[int]], + transform: np.ndarray, + method: str, +) -> None: + """Apply an existing transform to one moving tile and write to ``bw``. + Args: + br_mov: Moving ``BioReader``. + bw: Output ``BioWriter``. + tiles: ``(xm_bounds, ym_bounds, xt_bounds, yt_bounds)``. + shape: ``(tile_x, tile_y, x_crop_bounds, y_crop_bounds)``. + transform: Matrix from ``register_image``. + method: ``Projective`` or affine family. """ - - # Get the tile indices - Xm,Ym,Xt,Yt = tiles - + xm_bounds, ym_bounds, xt_bounds, yt_bounds = tiles + # Read the moving image tile - mov_tile = br_mov.read_image(X=[Xm[0],Xm[1]],Y=[Ym[0],Ym[1]],Z=[0,1],C=[0],T=[0]).squeeze() - + mov_tile = br_mov.read_image( + X=[xm_bounds[0], xm_bounds[1]], + Y=[ym_bounds[0], ym_bounds[1]], + Z=[0, 1], + C=[0], + T=[0], + ).squeeze() + # Get the image coordinates and shape - x,y,X_crop,Y_crop = shape - + tile_x, tile_y, x_crop_bounds, y_crop_bounds = shape + # Transform the moving image - if method=='Projective': - transformed_image = cv2.warpPerspective(mov_tile,transform,(Xt[1]-Xt[0],Yt[1]-Yt[0])) + if method == "Projective": + transformed_image = cv2.warpPerspective( + mov_tile, + transform, + (xt_bounds[1] - xt_bounds[0], yt_bounds[1] - yt_bounds[0]), + ) else: - transformed_image = cv2.warpAffine(mov_tile,transform,(Xt[1]-Xt[0],Yt[1]-Yt[0])) - + transformed_image = cv2.warpAffine( + mov_tile, + transform, + (xt_bounds[1] - xt_bounds[0], yt_bounds[1] - yt_bounds[0]), + ) + # Write the transformed image to the output file - bw.write_image(transformed_image[Y_crop[0]:Y_crop[1],X_crop[0]:X_crop[1],np.newaxis,np.newaxis,np.newaxis],[x],[y]) + bw.write_image( + transformed_image[ + y_crop_bounds[0] : y_crop_bounds[1], + x_crop_bounds[0] : x_crop_bounds[1], + np.newaxis, + np.newaxis, + np.newaxis, + ], + X=[tile_x], + Y=[tile_y], + ) + -if __name__=="__main__": - +if __name__ == "__main__": # Initialize the logger - logging.basicConfig(format='%(asctime)s - %(name)-8s - %(levelname)-8s - %(message)s', - datefmt='%d-%b-%y %H:%M:%S') + logging.basicConfig( + format="%(asctime)s - %(name)-8s - %(levelname)-8s - %(message)s", + datefmt="%d-%b-%y %H:%M:%S", + ) logger = logging.getLogger("image_registration.py") logger.setLevel(logging.INFO) # Setup the argument parsing logger.info("Parsing arguments...") - parser = argparse.ArgumentParser(prog='imageRegistration', description='This script registers an image collection') - parser.add_argument('--registrationString', dest='registration_string', type=str, required=True) - parser.add_argument('--similarTransformationString', dest='similar_transformation_string', type=str, required=True) - parser.add_argument('--outDir', dest='outDir', type=str, required=True) - parser.add_argument('--template', dest='template', type=str, required=True) - parser.add_argument('--method', dest='method', type=str, required=True) - - # parse the arguments + parser = argparse.ArgumentParser( + prog="imageRegistration", + description="This script registers an image collection", + ) + parser.add_argument( + "--registrationString", + dest="registration_string", + type=str, + required=True, + ) + parser.add_argument( + "--similarTransformationString", + dest="similar_transformation_string", + type=str, + required=True, + ) + parser.add_argument("--outDir", dest="outDir", type=str, required=True) + parser.add_argument("--template", dest="template", type=str, required=True) + parser.add_argument("--method", dest="method", type=str, required=True) + + # parse the arguments args = parser.parse_args() registration_string = args.registration_string similar_transformation_string = args.similar_transformation_string - outDir = args.outDir + output_dir = args.outDir template = args.template method = args.method - + # Set up the number of threads for each task - read_workers = max([cpu_count()//3,1]) - write_workers = max([cpu_count()-1,2]) - loop_workers = max([3*cpu_count()//4,2]) - + read_workers = max([cpu_count() // 3, 1]) + write_workers = max([cpu_count() - 1, 2]) + loop_workers = max([3 * cpu_count() // 4, 2]) + # extract filenames from registration_string and similar_transformation_string - registration_set=registration_string.split() - similar_transformation_set=similar_transformation_string.split() - - filename_len=len(template) - - # seperate the filename of the moving image from the complete path - moving_image_name=registration_set[1][-1*filename_len:] - + registration_set = registration_string.split() + similar_transformation_set = similar_transformation_string.split() + + filename_len = len(template) + # read and downscale reference image - br_ref = BioReader(registration_set[0],max_workers=write_workers) - scale_factor=get_scale_factor(br_ref.num_y(),br_ref.num_x()) - logger.info('Scale factor: {}'.format(scale_factor)) - - # intialize the scale factor and scale matrix(to be used to upscale the transformation matrices) - if method == 'Projective': - scale_matrix = np.array([[1,1,scale_factor],[1,1,scale_factor],[1/scale_factor,1/scale_factor,1]]) + br_ref = BioReader(registration_set[0], max_workers=write_workers) + scale_factor = get_scale_factor(br_ref.num_y(), br_ref.num_x()) + logger.info(f"Scale factor: {scale_factor}") + + # scale matrix used to upscale transformation matrices + if method == "Projective": + inv_sf = 1 / scale_factor + scale_matrix = np.array( + [ + [1, 1, scale_factor], + [1, 1, scale_factor], + [inv_sf, inv_sf, 1], + ], + ) else: - scale_matrix = np.array([[1/scale_factor,1/scale_factor,1],[1/scale_factor,1/scale_factor,1]]) - - logger.info('Reading and downscaling reference image: {}'.format(Path(registration_set[0]).name)) - reference_image_downscaled,max_val,min_val = get_scaled_down_images(br_ref,scale_factor,get_max=True) + inv_sf = 1 / scale_factor + row = [inv_sf, inv_sf, 1] + scale_matrix = np.array([row, row]) + + ref_name = Path(registration_set[0]).name + logger.info(f"Reading and downscaling reference image: {ref_name}") + ref_ds, max_val, min_val = get_scaled_down_images( + br_ref, + scale_factor, + get_max=True, + ) br_ref.max_workers = read_workers - + # read moving image - logger.info('Reading and downscaling moving image: {}'.format(Path(registration_set[1]).name)) - br_mov = BioReader(registration_set[1],max_workers=write_workers) - moving_image_downscaled = get_scaled_down_images(br_mov,scale_factor) + mov_name = Path(registration_set[1]).name + logger.info(f"Reading and downscaling moving image: {mov_name}") + br_mov = BioReader(registration_set[1], max_workers=write_workers) + moving_ds = get_scaled_down_images(br_mov, scale_factor) br_mov.max_workers = read_workers - + # calculate rough transformation between scaled down reference and moving image logger.info("calculating rough homography...") - Rough_Homography_Downscaled = get_transform(moving_image_downscaled, - reference_image_downscaled, - max_val, - min_val, - method) - + rough_homography_downscaled = get_transform( + moving_ds, + ref_ds, + max_val, + min_val, + method, + ) + # upscale the rough homography matrix logger.info("Inverting homography...") - if method=='Projective': - Rough_Homography_Upscaled=Rough_Homography_Downscaled*scale_matrix - homography_inverse=np.linalg.inv(Rough_Homography_Upscaled) + if method == "Projective": + rough_homography_upscaled = rough_homography_downscaled * scale_matrix + homography_inverse = np.linalg.inv(rough_homography_upscaled) else: - Rough_Homography_Upscaled=Rough_Homography_Downscaled - homography_inverse=cv2.invertAffineTransform(Rough_Homography_Downscaled) - + rough_homography_upscaled = rough_homography_downscaled + homography_inverse = cv2.invertAffineTransform(rough_homography_downscaled) + # Initialize the output file - bw = BioWriter(str(Path(outDir).joinpath(Path(registration_set[1]).name)),metadata=br_mov.read_metadata(),max_workers=write_workers) + out_path = str(Path(output_dir).joinpath(Path(registration_set[1]).name)) + bw = BioWriter( + out_path, + metadata=br_mov.read_metadata(), + max_workers=write_workers, + ) bw.num_x(br_ref.num_x()) bw.num_y(br_ref.num_y()) bw.num_z(1) bw.num_c(1) bw.num_t(1) - + # transformation variables reg_shape = [] reg_tiles = [] reg_homography = [] - + # Loop through image tiles and start threads logger.info("Starting threads...") threads = [] first_tile = True with ThreadPoolExecutor(max_workers=loop_workers) as executor: - for x in range(0,br_ref.num_x(),2048): - for y in range(0,br_ref.num_y(),2048): - + for x in range(0, br_ref.num_x(), 2048): + for y in range(0, br_ref.num_y(), 2048): # Get reference/template image coordinates - Xt = [np.max([0,x-1024]),np.min([br_ref.num_x(),x+2048+1024])] - Yt = [np.max([0,y-1024]),np.min([br_ref.num_y(),y+2048+1024])] - + nx, ny = br_ref.num_x(), br_ref.num_y() + xt_hi = min(nx, x + 2048 + 1024) + yt_hi = min(ny, y + 2048 + 1024) + xt_bounds = [max(0, x - 1024), xt_hi] + yt_bounds = [max(0, y - 1024), yt_hi] + # Use the rough homography to get coordinates in the moving image - coords = np.array([[Xt[0],Xt[0],Xt[1],Xt[1]], - [Yt[0],Yt[1],Yt[1],Yt[0]], - [1,1,1,1]], - dtype=np.float64) - - coords = np.matmul(homography_inverse,coords) - - mins = np.min(coords,axis=1) - maxs = np.max(coords,axis=1) - - Xm = [int(np.floor(np.max([mins[0],0]))), - int(np.ceil(np.min([maxs[0],br_mov.num_x()])))] - Ym = [int(np.floor(np.max([mins[1],0]))), - int(np.ceil(np.min([maxs[1],br_mov.num_y()])))] - - reg_tiles.append((Xm,Ym,Xt,Yt)) - + coords = np.array( + [ + [xt_bounds[0], xt_bounds[0], xt_bounds[1], xt_bounds[1]], + [yt_bounds[0], yt_bounds[1], yt_bounds[1], yt_bounds[0]], + [1, 1, 1, 1], + ], + dtype=np.float64, + ) + + coords = np.matmul(homography_inverse, coords) + + mins = np.min(coords, axis=1) + maxs = np.max(coords, axis=1) + + xm_bounds = [ + int(np.floor(np.max([mins[0], 0]))), + int(np.ceil(np.min([maxs[0], br_mov.num_x()]))), + ] + ym_bounds = [ + int(np.floor(np.max([mins[1], 0]))), + int(np.ceil(np.min([maxs[1], br_mov.num_y()]))), + ] + + reg_tiles.append((xm_bounds, ym_bounds, xt_bounds, yt_bounds)) + # Get cropping dimensions - X_crop = [1024 if Xt[0] > 0 else 0] - X_crop.append(2048+X_crop[0] if Xt[1]-Xt[0] >= 3072 else Xt[1]-Xt[0]+X_crop[0]) - Y_crop = [1024 if Yt[0] > 0 else 0] - Y_crop.append(2048+Y_crop[0] if Yt[1]-Yt[0] >= 3072 else Yt[1]-Yt[0]+Y_crop[0]) - reg_shape.append((x,y,X_crop,Y_crop)) - + x_crop_bounds = [1024 if xt_bounds[0] > 0 else 0] + x_crop_bounds.append( + 2048 + x_crop_bounds[0] + if xt_bounds[1] - xt_bounds[0] >= TILE_MERGE_THRESHOLD + else xt_bounds[1] - xt_bounds[0] + x_crop_bounds[0], + ) + y_crop_bounds = [1024 if yt_bounds[0] > 0 else 0] + y_crop_bounds.append( + 2048 + y_crop_bounds[0] + if yt_bounds[1] - yt_bounds[0] >= TILE_MERGE_THRESHOLD + else yt_bounds[1] - yt_bounds[0] + y_crop_bounds[0], + ) + reg_shape.append((x, y, x_crop_bounds, y_crop_bounds)) + # Start a thread to register the tiles - threads.append(executor.submit(register_image,br_ref,br_mov,bw,Xt,Yt,Xm,Ym,x,y,X_crop,Y_crop,max_val,min_val,method)) - + threads.append( + executor.submit( + register_image, + br_ref, + br_mov, + bw, + xt_bounds, + yt_bounds, + xm_bounds, + ym_bounds, + x, + y, + x_crop_bounds, + y_crop_bounds, + max_val, + min_val, + method, + rough_homography_upscaled, + ), + ) + # Bioformats require the first tile be written before any other tile if first_tile: - logger.info('Waiting for first_tile to finish...') + logger.info("Waiting for first_tile to finish...") first_tile = False threads[0].result() - + # Wait for threads to finish, track progress for thread_num in range(len(threads)): if thread_num % 10 == 0: - logger.info('Registration progress: {:6.2f}%'.format(100*thread_num/len(threads))) + pct = 100 * thread_num / len(threads) + logger.info(f"Registration progress: {pct:6.2f}%") reg_homography.append(threads[thread_num].result()) - + # Close the image bw.close_image() - logger.info('Registration progress: {:6.2f}%'.format(100.0)) - - # iterate across all images which have the similar transformation as the moving image above + logger.info(f"Registration progress: {100.0:6.2f}%") + + # remaining images sharing the moving image's transform for moving_image_path in similar_transformation_set: - # seperate image name from the path to it - moving_image_name=moving_image_path[-1*filename_len:] - - logger.info('Applying registration to image: {}'.format(moving_image_name)) - - br_mov = BioReader(moving_image_path,max_workers=read_workers) - - bw = BioWriter(str(Path(outDir).joinpath(moving_image_name)), metadata=br_mov.read_metadata(),max_workers=write_workers) + moving_image_name = moving_image_path[-1 * filename_len :] + + logger.info(f"Applying registration to image: {moving_image_name}") + + br_mov = BioReader(moving_image_path, max_workers=read_workers) + + bw = BioWriter( + str(Path(output_dir).joinpath(moving_image_name)), + metadata=br_mov.read_metadata(), + max_workers=write_workers, + ) bw.num_x(br_ref.num_x()) bw.num_y(br_ref.num_y()) bw.num_z(1) bw.num_c(1) bw.num_t(1) - + # Apply transformation to remaining images - logger.info('Transformation progress: {:5.2f}%'.format(0.0)) + logger.info(f"Transformation progress: {0.0:5.2f}%") threads = [] - with ThreadPoolExecutor(loop_workers) as executor: + with ThreadPoolExecutor(max_workers=loop_workers) as executor: first_tile = True - for tile,shape,transform in zip(reg_tiles,reg_shape,reg_homography): - + for tile, shape, transform in zip(reg_tiles, reg_shape, reg_homography): # Start transformation threads - threads.append(executor.submit(apply_transform,br_mov,bw,tile,shape,transform,method)) - + threads.append( + executor.submit( + apply_transform, + br_mov, + bw, + tile, + shape, + transform, + method, + ), + ) + # The first tile must be written before all other tiles if first_tile: first_tile = False threads[0].result() - # Wait for threads to finish and track progress + # Wait for threads to finish and track progress for thread_num in range(len(threads)): if thread_num % 10 == 0: - logger.info('Transformation progress: {:6.2f}%'.format(100*thread_num/len(threads))) + tpct = 100 * thread_num / len(threads) + logger.info(f"Transformation progress: {tpct:6.2f}%") threads[thread_num].result() - logger.info('Transformation progress: {:6.2f}%'.format(100.0)) - + logger.info(f"Transformation progress: {100.0:6.2f}%") + bw.close_image() - \ No newline at end of file diff --git a/transforms/images/polus-image-registration-plugin/src/log4j.properties b/transforms/images/polus-image-registration-plugin/src/log4j.properties index b6682b3d8..bdcc5504b 100644 --- a/transforms/images/polus-image-registration-plugin/src/log4j.properties +++ b/transforms/images/polus-image-registration-plugin/src/log4j.properties @@ -9,4 +9,4 @@ log4j.appender.A1=org.apache.log4j.ConsoleAppender # A1 uses PatternLayout. log4j.appender.A1.layout=org.apache.log4j.PatternLayout -log4j.appender.A1.layout.ConversionPattern=%-4r [%t] %-5p %c %x - %m%n \ No newline at end of file +log4j.appender.A1.layout.ConversionPattern=%-4r [%t] %-5p %c %x - %m%n diff --git a/transforms/images/polus-image-registration-plugin/src/main.py b/transforms/images/polus-image-registration-plugin/src/main.py index 8c6e05a92..48a596b31 100644 --- a/transforms/images/polus-image-registration-plugin/src/main.py +++ b/transforms/images/polus-image-registration-plugin/src/main.py @@ -1,78 +1,151 @@ -import argparse, logging, subprocess -import numpy as np -from pathlib import Path -from parser import parse_collection -import shutil +"""CLI entry: parse collections and delegate to ``image_registration`` workers.""" +from __future__ import annotations +import argparse +import logging +import shutil +import subprocess +import sys +from parser import parse_collection +from pathlib import Path -if __name__=="__main__": +if __name__ == "__main__": # Initialize the logger - logging.basicConfig(format='%(asctime)s - %(name)-8s - %(levelname)-8s - %(message)s', - datefmt='%d-%b-%y %H:%M:%S') + logging.basicConfig( + format="%(asctime)s - %(name)-8s - %(levelname)-8s - %(message)s", + datefmt="%d-%b-%y %H:%M:%S", + ) logger = logging.getLogger("main") logger.setLevel(logging.INFO) # Setup the argument parsing logger.info("Parsing arguments...") - parser = argparse.ArgumentParser(prog='main', description='This plugin registers an image collection') - parser.add_argument('--filePattern', dest='filePattern', type=str, help='Filename pattern used to separate data', required=True) - parser.add_argument('--inpDir', dest='inpDir', type=str,help='Input image collection to be processed by this plugin', required=True) - parser.add_argument('--registrationVariable', dest='registrationVariable', type=str, help='variable to help identify which images need to be registered to each other', required=True) - parser.add_argument('--template', dest='template', type=str, help='Template image to be used for image registration', required=True) - parser.add_argument('--TransformationVariable', dest='TransformationVariable', type=str,help='variable to help identify which images have similar transformation', required=True) - parser.add_argument('--outDir', dest='outDir', type=str, help='Output collection', required=True) - parser.add_argument('--method', dest='method', type=str, help='projective, affine, or partialaffine', required=True) - - # Parse the arguments + parser = argparse.ArgumentParser( + prog="main", + description="This plugin registers an image collection", + ) + parser.add_argument( + "--filePattern", + dest="filePattern", + type=str, + help="Filename pattern used to separate data", + required=True, + ) + parser.add_argument( + "--inpDir", + dest="inpDir", + type=str, + help="Input image collection to be processed by this plugin", + required=True, + ) + parser.add_argument( + "--registrationVariable", + dest="registrationVariable", + type=str, + help=( + "variable to help identify which images need to be registered to each other" + ), + required=True, + ) + parser.add_argument( + "--template", + dest="template", + type=str, + help="Template image to be used for image registration", + required=True, + ) + parser.add_argument( + "--TransformationVariable", + dest="TransformationVariable", + type=str, + help="variable to help identify which images have similar transformation", + required=True, + ) + parser.add_argument( + "--outDir", + dest="outDir", + type=str, + help="Output collection", + required=True, + ) + parser.add_argument( + "--method", + dest="method", + type=str, + help="projective, affine, or partialaffine", + required=True, + ) + + # Parse the arguments args = parser.parse_args() - filePattern = args.filePattern - logger.info('filePattern = {}'.format(filePattern)) - - inpDir = args.inpDir + file_pattern = args.filePattern + logger.info(f"filePattern = {file_pattern}") + + input_dir = args.inpDir # check if images folder is present in the input directory - if (Path.is_dir(Path(inpDir).joinpath('images'))): - inpDir= str(Path(inpDir).joinpath('images')) - - logger.info('inpDir = {}'.format(inpDir)) - registrationVariable = args.registrationVariable - logger.info('registrationVariable = {}'.format(registrationVariable)) + if Path.is_dir(Path(input_dir).joinpath("images")): + input_dir = str(Path(input_dir).joinpath("images")) + + logger.info(f"inpDir = {input_dir}") + registration_variable = args.registrationVariable + logger.info(f"registrationVariable = {registration_variable}") template = args.template - logger.info('template = {}'.format(template)) - TransformationVariable = args.TransformationVariable - logger.info('TransformationVariable = {}'.format(TransformationVariable)) - outDir = args.outDir - logger.info('outDir = {}'.format(outDir)) + logger.info(f"template = {template}") + transformation_variable = args.TransformationVariable + logger.info(f"TransformationVariable = {transformation_variable}") + output_dir = args.outDir + logger.info(f"outDir = {output_dir}") method = args.method - logger.info('method = {}'.format(method)) - - # get template image path - template_image_path=str(Path(inpDir).joinpath(template).absolute()) - + logger.info(f"method = {method}") + + # get template image path + template_image_path = str(Path(input_dir).joinpath(template).absolute()) + # filename len - filename_len= len(template) - + filename_len = len(template) + # parse the input collection - logger.info('Parsing the input collection and getting registration_dictionary') - registration_dictionary=parse_collection(inpDir,filePattern,registrationVariable, TransformationVariable, template_image_path) - - logger.info('Iterating over registration_dictionary....') - for registration_set,similar_transformation_set in registration_dictionary.items(): - + logger.info("Parsing the input collection and getting registration_dictionary") + registration_dictionary = parse_collection( + input_dir, + file_pattern, + registration_variable, + transformation_variable, + template_image_path, + ) + + logger.info("Iterating over registration_dictionary....") + script_path = Path(__file__).resolve().parent / "image_registration.py" + for registration_set, similar_transformation_set in registration_dictionary.items(): # registration_dictionary consists of set of already registered images as well - if registration_set[0]==registration_set[1]: - similar_transformation_set=similar_transformation_set.tolist() - similar_transformation_set.append(registration_set[0]) - for image_path in similar_transformation_set: - image_name=image_path[-1*filename_len:] - logger.info('Copying image {} to output directory'.format(image_name)) - shutil.copy2(image_path,str(Path(outDir).joinpath(image_name).absolute())) + if registration_set[0] == registration_set[1]: + paths_to_copy = list(similar_transformation_set) + paths_to_copy.append(registration_set[0]) + for image_path in paths_to_copy: + image_name = image_path[-1 * filename_len :] + logger.info(f"Copying image {image_name} to output directory") + dest = Path(output_dir).joinpath(image_name).absolute() + shutil.copy2(image_path, str(dest)) continue - + # concatenate lists into a string to pass as an argument to argparse - registration_string=' '.join(registration_set) - similar_transformation_string=' '.join(similar_transformation_set) + registration_string = " ".join(registration_set) + similar_transformation_string = " ".join(similar_transformation_set) - # open subprocess image_registration.py - registration = subprocess.Popen("python3 image_registration.py --registrationString '{}' --similarTransformationString '{}' --outDir '{}' --template '{}' --method '{}'".format(registration_string,similar_transformation_string,outDir,template,method), shell=True ) - registration.wait() - \ No newline at end of file + subprocess.run( + [ # noqa: S603 + sys.executable, + str(script_path), + "--registrationString", + registration_string, + "--similarTransformationString", + similar_transformation_string, + "--outDir", + output_dir, + "--template", + template, + "--method", + method, + ], + check=True, + ) diff --git a/transforms/images/polus-image-registration-plugin/src/parser.py b/transforms/images/polus-image-registration-plugin/src/parser.py index 5b7560fa5..1057b4a2b 100755 --- a/transforms/images/polus-image-registration-plugin/src/parser.py +++ b/transforms/images/polus-image-registration-plugin/src/parser.py @@ -1,117 +1,108 @@ -from filepattern import FilePattern, parse_directory,get_matching,parse_filename,get_regex -import os +"""Build registration groups from a filepattern-indexed image directory.""" +from __future__ import annotations + import itertools + import numpy as np +from filepattern import FilePattern +from filepattern import get_regex +from filepattern import parse_directory -def parse_collection(directory_path,file_pattern,registration_variable, similarity_variable, template_image): - - """ - This function parses the input directory and returns a dictionary. Each key in the dictionary is a tuple - consisting of a template and a moving image. The value corresponding to each key is a list of images that have - similar transformation as the moving image. - Note: The code produces the expected output when len(registration_variable)==len(similarity_variable)==1. - The code will NOT spit out an error when the more than one variable is passed as registration or - similarity variable, but additional testing needs to be done to validate the script for this usecase. - - inputs : - directory_path: path to the input collection - file_pattern: file name pattern of the input images - registration_variable : variable to help determine the set of moving and template images - similarity variable: variable to help determine the set of images having a similar transformation - corresponding to each set of moving and template images - template_image: name of a template image - - outputs : result_dic - example of result_dic is shown below - - result_dic = { (template_img1 : moving_img1) : [set1_img1,set1_img2, set1_img3....], - (template_img2 : moving_img2) : [set2_img1,set2_img2, set2_img3....], - . - . } - - - +def parse_collection( + directory_path: str, + file_pattern: str, + registration_variable: str, + similarity_variable: str, + _template_image: str, +) -> dict[tuple[str, str], np.ndarray]: + r"""Parse the input directory into a registration / similar-transform map. + + Each key is ``(template_path, moving_path)``; the value lists paths that share + the moving image's transform. + + Note: + Tested with ``len(registration_variable) == len(similarity_variable) == 1``. + Other lengths are not validated. + + Args: + directory_path: Path to the input collection. + file_pattern: Filename pattern of the input images. + registration_variable: Variable that groups template vs moving images. + similarity_variable: Variable for images with the same transform as moving. + _template_image: Reserved for future use (template basename). + + Returns: + ``result_dic`` mapping ``(template, moving)`` to an array of related paths. """ - - # Predefined variables order - #var_order = 'rtczyx' - # get all variables in the file pattern - _,variables=get_regex(file_pattern) - + _, variables = get_regex(file_pattern) + # get variables except the registration and similarity variable - moving_variables=[var for var in variables if var not in registration_variable and var not in similarity_variable] - + moving_variables = [ + var + for var in variables + if var not in registration_variable and var not in similarity_variable + ] + # uvals is dictionary with all the possible variables as key - # corresponding to each key is a list of all values which that variable can take for the input collection - _, uvals=parse_directory(directory_path,file_pattern) - - parser_object=FilePattern(directory_path,file_pattern) - - image_set=[] - - # extract the index values from uvals for each variable in moving_variables - moving_variables_set=[uvals[var] for var in moving_variables] - - # iterate over the similar transformation variables - # Code produced expected output when len(registration_variable)==len(similarity_variable)==1 - # refer to function description + # corresponding to each key is a list of all values which that variable can take + _, uvals = parse_directory(directory_path, file_pattern) + + parser_object = FilePattern(directory_path, file_pattern) + + image_set = [] + + # extract the index values from uvals for each variable in moving_variables + moving_variables_set = [uvals[var] for var in moving_variables] + + # iterate over the similar transformation variables; + # expected output when len(registration_variable)==len(similarity_variable)==1 for char in similarity_variable: - # append the variable to the moving variable set - moving_variables.append(char) - + moving_variables.append(char) + # iterate over all possible index values of the similar transf. variable for ind in uvals[char]: - registration_set=[] - + registration_set = [] + # append the fixed value of the index to the moving variables set - moving_variables_set.append([ind]) - - # get all the possible combinations of the index values in the moving variables set - registration_indices_combinations=list(itertools.product(*moving_variables_set)) - all_dicts=[] - - # iterate over all combinations and create a dictionary for each combination - # the dictionary is of the form {'C'=1, 'X'=2...etc} which can be used as an input - # to the get_matching() function - for index_comb in registration_indices_combinations: - inter_dict={} + moving_variables_set.append([ind]) + + # all combinations of index values in the moving variables set + reg_idx_combos = list(itertools.product(*moving_variables_set)) + all_dicts = [] + + # dicts like {'C': 1, 'X': 2, ...} for get_matching() + for index_comb in reg_idx_combos: + inter_dict = {} for i in range(len(moving_variables)): - inter_dict.update({moving_variables[i].upper():index_comb[i]}) + inter_dict.update({moving_variables[i].upper(): index_comb[i]}) # store all dictionaries - all_dicts.append(inter_dict) - + all_dicts.append(inter_dict) + # iterate over all dictionaries for reg_dict in all_dicts: - intermediate_set=[] - # use get_matching function to get all filenames with defined variable values in the dictionary - files=parser_object.get_matching(**reg_dict) - - # files is a list of dictionaries + intermediate_set = [] + files = parser_object.get_matching(**reg_dict) + + # files is a list of dictionaries for file_dict in files: - intermediate_set.append(file_dict['file']) - registration_set.append(intermediate_set) - - # delete the fixed index value of the similar transf. variable to prepare for the next iteration + intermediate_set.append(file_dict["file"]) + registration_set.append(intermediate_set) + + # drop fixed similar-transform index before next iteration moving_variables_set.pop(-1) image_set.append(registration_set) - - - # parse image set to form the result dictionary - result_dic={} - old_set=np.array(image_set) + + # parse image set to form the result dictionary + result_dic = {} + old_set = np.array(image_set) for j in range(old_set.shape[1]): - inter=old_set[:,j,:] + inter = old_set[:, j, :] for k in range(inter.shape[1]): - ky=(inter[0,0],inter[0,k]) - items=inter[1:,k] - result_dic.update({ky:items}) - - return result_dic - + ky = (inter[0, 0], inter[0, k]) + items = inter[1:, k] + result_dic.update({ky: items}) - - - \ No newline at end of file + return result_dic