From 1c30ea4a8f1d6695860b5bff666e99c38648fc9f Mon Sep 17 00:00:00 2001
From: James Ball <37094972+PatBall1@users.noreply.github.com>
Date: Tue, 17 Sep 2024 12:24:43 +0100
Subject: [PATCH] Jb/july24 (#148)
* parallelise tiling
* multispectral compatibility
* updated multi-class functionality
---
.gitattributes | 1 -
README.md | 8 +-
detectree2/data_loading/custom.py | 74 +++
detectree2/models/train.py | 756 ++++++++++++++++++-------
detectree2/preprocessing/tiling.py | 639 ++++++++++++---------
detectree2/tests/test_preprocessing.py | 4 +-
docs/source/index.rst | 1 +
docs/source/tutorial.rst | 432 +++++++++++++-
docs/source/tutorial_multi.rst | 166 +++++-
report/figures/IoU_AP.png | Bin 0 -> 251986 bytes
report/figures/train_val_loss.png | Bin 0 -> 41469 bytes
report/figures/val_AP50.png | Bin 0 -> 22751 bytes
setup.py | 4 +-
13 files changed, 1562 insertions(+), 523 deletions(-)
create mode 100644 detectree2/data_loading/custom.py
create mode 100644 report/figures/IoU_AP.png
create mode 100644 report/figures/train_val_loss.png
create mode 100644 report/figures/val_AP50.png
diff --git a/.gitattributes b/.gitattributes
index ec4a626f..e69de29b 100644
--- a/.gitattributes
+++ b/.gitattributes
@@ -1 +0,0 @@
-*.pth filter=lfs diff=lfs merge=lfs -text
diff --git a/README.md b/README.md
index 2f07f770..0958d783 100644
--- a/README.md
+++ b/README.md
@@ -10,18 +10,14 @@
-Python package for automatic tree crown delineation based on Mask R-CNN. Pre-trained models can be picked in the [`model_garden`](https://github.com/PatBall1/detectree2/tree/master/model_garden).
-A tutorial on how to prepare data, train models and make predictions is available [here](https://patball1.github.io/detectree2/tutorial.html). For questions, collaboration proposals and requests for data email [James Ball](mailto:ball.jgc@gmail.com). Some example data is available for download [here](https://doi.org/10.5281/zenodo.8136161).
+Python package for automatic tree crown delineation in aerial RGB and multispectral imagery based on Mask R-CNN. Pre-trained models can be picked in the [`model_garden`](https://github.com/PatBall1/detectree2/tree/master/model_garden).
+A tutorial on how to prepare data, train models and make predictions is available [here](https://patball1.github.io/detectree2/tutorial.html). For questions, collaboration proposals and requests for data email [James Ball](mailto:ball.jgc@gmail.com). Some example data is available to download [here](https://doi.org/10.5281/zenodo.8136161).
Detectree2是一个基于Mask R-CNN的自动树冠检测与分割的Python包。您可以在[`model_garden`](https://github.com/PatBall1/detectree2/tree/master/model_garden)中选择预训练模型。[这里](https://patball1.github.io/detectree2/tutorial.html)提供了如何准备数据、训练模型和进行预测的教程。如果有任何问题，合作提案或者需要样例数据，可以邮件联系[James Ball](mailto:ball.jgc@gmail.com)。一些示例数据可以在[这里](https://doi.org/10.5281/zenodo.8136161)下载。
| 
| Code developed by James Ball, Seb Hickman, Thomas Koay, Oscar Jiang, Luran Wang, Panagiotis Ioannou, James Hinton and Matthew Archer in the [Forest Ecology and Conservation Group](https://coomeslab.org/) at the University of Cambridge. The Forest Ecology and Conservation Group is led by Professor David Coomes and is part of the University of Cambridge [Conservation Research Institute](https://www.conservation.cam.ac.uk/). |
| :---: | :--- |
-
-> [!NOTE]
-> To save bandwidth, trained models have been moved to [Zenodo](https://zenodo.org/records/10522461). Download models directly with `wget` or equivalent.
-
## Citation
diff --git a/detectree2/data_loading/custom.py b/detectree2/data_loading/custom.py
new file mode 100644
index 00000000..81a1becb
--- /dev/null
+++ b/detectree2/data_loading/custom.py
@@ -0,0 +1,74 @@
+import cv2
+import detectron2.data.transforms as T
+import numpy as np
+import rasterio
+import torch
+from detectron2.structures import BitMasks, BoxMode, Instances
+from torch.utils.data import Dataset
+
+
+class CustomTIFFDataset(Dataset):
+ def __init__(self, annotations, transforms=None):
+ """
+ Args:
+ annotations (list): List of dictionaries containing image file paths and annotations.
+ transforms (callable, optional): Optional transform to be applied on a sample.
+ """
+ self.annotations = annotations
+ self.transforms = transforms
+
+ def __len__(self):
+ return len(self.annotations)
+
+ def __getitem__(self, idx):
+ # Load the TIFF image
+ img_info = self.annotations[idx]
+ with rasterio.open(img_info['file_name']) as src:
+ # Read all bands (assuming they are all needed)
+ image = src.read()
+ # Normalize or rescale if necessary
+ image = image.astype(np.float32) / 255.0 # Example normalization
+ # If the number of bands is not 3, reduce to 3 or handle accordingly
+ #if image.shape[0] > 3:
+ # image = image[:3, :, :] # Taking the first 3 bands (e.g., RGB)
+ # Convert to HWC format expected by Detectron2
+ #image = np.transpose(image, (1, 2, 0))
+
+ # Prepare annotations (this part needs to be adapted to your specific annotations)
+ target = {
+ "image_id": idx,
+ "annotations": img_info["annotations"],
+ "width": img_info["width"],
+ "height": img_info["height"],
+ }
+
+ if self.transforms is not None:
+ augmentations = T.AugmentationList(self.transforms)
+ image, target = augmentations(image, target)
+
+ # Convert to Detectron2-compatible format
+ image = torch.as_tensor(image.astype("float32").transpose(2, 0, 1))
+ instances = self.get_detectron_instances(target)
+
+ return image, instances
+
+ def get_detectron_instances(self, target):
+ """
+ Converts annotations into Detectron2's format.
+ This example assumes annotations are in COCO format, and you'll need to adapt it for your needs.
+ """
+ boxes = [obj["bbox"] for obj in target["annotations"]]
+ boxes = torch.as_tensor(boxes, dtype=torch.float32)
+ boxes = BoxMode.convert(boxes, BoxMode.XYWH_ABS, BoxMode.XYXY_ABS)
+
+ # Create BitMasks from the binary mask data (assuming the mask is a binary numpy array)
+ masks = [obj["segmentation"] for obj in target["annotations"]] # Replace with actual mask loading
+ masks = BitMasks(torch.stack([torch.from_numpy(mask) for mask in masks]))
+
+ instances = Instances(
+ image_size=(target["height"], target["width"]),
+ gt_boxes=boxes,
+ gt_classes=torch.tensor([obj["category_id"] for obj in target["annotations"]], dtype=torch.int64),
+ gt_masks=masks
+ )
+ return instances
diff --git a/detectree2/models/train.py b/detectree2/models/train.py
index d7606793..5e404e4d 100644
--- a/detectree2/models/train.py
+++ b/detectree2/models/train.py
@@ -9,15 +9,18 @@
import logging
import os
import random
+import re
import time
from pathlib import Path
-from typing import Any, Dict, List
+from typing import Any, Dict, List, Optional
import cv2
import detectron2.data.transforms as T # noqa:N812
import detectron2.utils.comm as comm
import numpy as np
+import rasterio
import torch
+import torch.nn as nn
from detectron2 import model_zoo
from detectron2.checkpoint import DetectionCheckpointer # noqa:F401
from detectron2.config import get_cfg
@@ -28,6 +31,7 @@
build_detection_test_loader,
build_detection_train_loader,
)
+from detectron2.data import detection_utils as utils
from detectron2.engine import DefaultTrainer
from detectron2.engine.hooks import HookBase
from detectron2.evaluation import COCOEvaluator, verify_results
@@ -38,29 +42,144 @@
from detectron2.utils.logger import log_every_n_seconds
from detectron2.utils.visualizer import ColorMode, Visualizer
-# from IPython.display import display
-# from PIL import Image
+from detectree2.preprocessing.tiling import load_class_mapping
-class LossEvalHook(HookBase):
- """Do inference and get the loss metric.
+class FlexibleDatasetMapper(DatasetMapper):
+ """
+ A flexible dataset mapper that extends the standard DatasetMapper to handle
+ multi-band images and custom augmentations.
- Class to:
- - Do inference of dataset like an Evaluator does
- - Get the loss metric like the trainer does
- https://github.com/facebookresearch/detectron2/blob/master/detectron2/evaluation/evaluator.py
- https://github.com/facebookresearch/detectron2/blob/master/detectron2/engine/train_loop.py
- See https://gist.github.com/ortegatron/c0dad15e49c2b74de8bb09a5615d9f6b
+ This class is designed to work with datasets that may contain images with
+ more than three channels (e.g., multispectral images) and allows for custom
+ augmentations to be applied. It also handles semantic segmentation data if
+ provided in the dataset.
+
+ Args:
+ cfg (CfgNode): Configuration object containing dataset and model configurations.
+ is_train (bool): Flag indicating whether the mapper is being used for training. Default is True.
+ augmentations (list, optional): List of augmentations to be applied. Default is an empty list.
Attributes:
- model: model to train
- period: number of iterations between evaluations
- data_loader: data loader to use for evaluation
- patience: number of evaluation periods to wait for improvement
+ cfg (CfgNode): Stores the configuration object for later use.
+ is_train (bool): Indicates whether the mapper is in training mode.
+ logger (Logger): Logger instance for logging messages.
+ """
+ def __init__(self, cfg, is_train=True, augmentations=None):
+ if augmentations is None:
+ augmentations = []
+
+ # Initialize the base DatasetMapper class with provided parameters
+ super().__init__(
+ is_train=is_train,
+ augmentations=augmentations,
+ image_format=cfg.INPUT.FORMAT,
+ use_instance_mask=cfg.MODEL.MASK_ON,
+ use_keypoint=cfg.MODEL.KEYPOINT_ON,
+ instance_mask_format=cfg.INPUT.MASK_FORMAT,
+ keypoint_hflip_indices=None,
+ precomputed_proposal_topk=None,
+ recompute_boxes=False
+ )
+ self.cfg = cfg
+ self.is_train = is_train
+ self.logger = logging.getLogger(__name__)
+ mode = "training" if is_train else "inference"
+ self.logger.info(f"[FlexibleDatasetMapper] Augmentations used in {mode}: {augmentations}")
+
+ def __call__(self, dataset_dict):
+ """
+ Process a single dataset dictionary, applying the necessary transformations and augmentations.
+
+ Args:
+ dataset_dict (dict): A dictionary containing data for a single dataset item, including
+ file names and metadata.
+
+ Returns:
+ dict: The processed dataset dictionary, or None if there was an error.
+ """
+ if dataset_dict is None:
+ self.logger.warning("Received None for dataset_dict, skipping this entry.")
+ return None
+
+ if self.cfg.IMGMODE == "rgb":
+ return super().__call__(dataset_dict)
+
+ try:
+ # Handle multi-band image loading using rasterio
+ with rasterio.open(dataset_dict["file_name"]) as src:
+ img = src.read()
+ if img is None:
+ raise ValueError(f"Image data is None for file: {dataset_dict['file_name']}")
+ # Transpose image dimensions to match expected format (H, W, C)
+ img = np.transpose(img, (1, 2, 0)).astype("float32")
+
+ # Size check similar to utils.check_image_size
+ if img.shape[:2] != (dataset_dict.get("height"), dataset_dict.get("width")):
+ self.logger.warning(
+ f"""Image size {img.shape[:2]} does not match expected size {(dataset_dict.get('height'),
+ dataset_dict.get('width'))}."""
+ )
+
+ # Otherwise, handle custom multi-band logic
+ aug_input = T.AugInput(img)
+ transforms = self.augmentations(aug_input) # Apply the augmentations
+ img = aug_input.image
+
+ dataset_dict["image"] = torch.as_tensor(np.ascontiguousarray(img.transpose(2, 0, 1)))
+
+ # Handle semantic segmentation if present
+ if "sem_seg_file_name" in dataset_dict:
+ sem_seg_gt = utils.read_image(dataset_dict.pop("sem_seg_file_name"), "L").squeeze(2)
+ dataset_dict["sem_seg"] = torch.as_tensor(sem_seg_gt.astype("long"))
+
+ if not self.is_train:
+ # If not in training mode, remove annotations and segmentation file names
+ dataset_dict.pop("annotations", None)
+ dataset_dict.pop("sem_seg_file_name", None)
+ return dataset_dict
+
+ if "annotations" in dataset_dict:
+ # Apply the transformations to the annotations
+ self._transform_annotations(dataset_dict, transforms, img.shape[:2])
+
+ return dataset_dict
+
+ except Exception as e:
+ file_name = dataset_dict.get('file_name', 'unknown') if dataset_dict else 'unknown'
+ self.logger.error(f"Error processing {file_name}: {e}")
+ return None
+
+
+class LossEvalHook(HookBase):
"""
+ A custom hook for evaluating loss during training and managing model checkpoints based on evaluation metrics.
+
+ This hook is designed to:
+ - Perform inference on a dataset similarly to an Evaluator.
+ - Calculate and log the loss metric during training.
+ - Save the best model checkpoint based on a specified evaluation metric (e.g., AP50).
+ - Implement early stopping if the evaluation metric does not improve over a specified number of evaluations.
+ Attributes:
+ _model: The model to evaluate.
+ _period: Number of iterations between evaluations.
+ _data_loader: The data loader used for evaluation.
+ patience: Number of evaluation periods to wait before early stopping.
+ iter: Tracks the number of evaluations since the last improvement in the evaluation metric.
+ max_ap: The best evaluation metric (e.g., AP50) achieved during training.
+ best_iter: The iteration at which the best evaluation metric was achieved.
+ """
def __init__(self, eval_period, model, data_loader, patience):
- """Inits LossEvalHook."""
+ """
+ Initialize the LossEvalHook.
+
+ Args:
+ eval_period (int): The number of iterations between evaluations.
+ model (torch.nn.Module): The model to evaluate.
+ data_loader (torch.utils.data.DataLoader): The data loader for evaluation.
+ patience (int): The number of evaluation periods to wait for improvement before early stopping.
+ """
self._model = model
self._period = eval_period
self._data_loader = data_loader
@@ -70,10 +189,14 @@ def __init__(self, eval_period, model, data_loader, patience):
self.best_iter = 0
def _do_loss_eval(self):
- """Copying inference_on_dataset from evaluator.py.
+ """
+ Perform inference on the dataset and calculate the average loss.
+
+ This method is adapted from `inference_on_dataset` in Detectron2's evaluator.
+ It also calculates and logs the AP50 metric and updates the best model checkpoint if needed.
Returns:
- _type_: _description_
+ list: A list of loss values for each batch in the dataset.
"""
total = len(self._data_loader)
num_warmup = min(5, total - 1)
@@ -83,6 +206,7 @@ def _do_loss_eval(self):
losses = []
for idx, inputs in enumerate(self._data_loader):
if idx == num_warmup:
+ # Reset the start time after the warm-up phase
start_time = time.perf_counter()
total_compute_time = 0
start_compute_time = time.perf_counter()
@@ -92,6 +216,7 @@ def _do_loss_eval(self):
iters_after_start = idx + 1 - num_warmup * int(idx >= num_warmup)
seconds_per_img = total_compute_time / iters_after_start
if idx >= num_warmup * 2 or seconds_per_img > 5:
+ # Log progress and estimated time remaining
total_seconds_per_img = (time.perf_counter() - start_time) / iters_after_start
eta = datetime.timedelta(seconds=int(total_seconds_per_img * (total - idx - 1)))
log_every_n_seconds(
@@ -100,11 +225,13 @@ def _do_loss_eval(self):
str(eta)),
n=5,
)
+ # Calculate loss for the current batch
loss_batch = self._get_loss(inputs)
losses.append(loss_batch)
+
mean_loss = np.mean(losses)
- # print(self.trainer.cfg.DATASETS.TEST)
- # Combine the AP50s of the different datasets
+
+ # Calculate the average AP50 across datasets if multiple datasets are used for testing
if len(self.trainer.cfg.DATASETS.TEST) > 1:
APs = []
for dataset in self.trainer.cfg.DATASETS.TEST:
@@ -112,7 +239,10 @@ def _do_loss_eval(self):
AP = sum(APs) / len(APs)
else:
AP = self.trainer.test(self.trainer.cfg, self.trainer.model)["segm"]["AP50"]
- print("Av. AP50 =", AP)
+
+ print("Av. segm AP50 =", AP)
+
+ # Store the calculated loss and AP50 in the trainer's storage
self.trainer.APs.append(AP)
self.trainer.storage.put_scalar("validation_loss", mean_loss)
self.trainer.storage.put_scalar("validation_ap", AP)
@@ -121,15 +251,17 @@ def _do_loss_eval(self):
return losses
def _get_loss(self, data):
- """Calculate loss in train_loop.
+ """
+ Compute the loss for a given batch of data.
Args:
- data (_type_): _description_
+ data (dict): A batch of input data.
Returns:
- _type_: _description_
+ float: The total loss for the batch.
"""
metrics_dict = self._model(data)
+ # Detach and move to CPU for logging
metrics_dict = {
k: v.detach().cpu().item() if isinstance(v, torch.Tensor) else float(v)
for k, v in metrics_dict.items()
@@ -138,57 +270,74 @@ def _get_loss(self, data):
return total_losses_reduced
def after_step(self):
+ """
+ Hook to be called after each training iteration to evaluate the model and manage checkpoints.
+
+ - Evaluates the model at regular intervals.
+ - Saves the best model checkpoint based on the AP50 metric.
+ - Implements early stopping if the AP50 does not improve after a set number of evaluations.
+ """
next_iter = self.trainer.iter + 1
is_final = next_iter == self.trainer.max_iter
if is_final or (self._period > 0 and next_iter % self._period == 0):
self._do_loss_eval()
+ # Check if the current AP50 is the best so far
if self.max_ap < self.trainer.APs[-1]:
self.iter = 0
self.max_ap = self.trainer.APs[-1]
+ # Save the current best model
self.trainer.checkpointer.save("model_" + str(len(self.trainer.APs)))
self.best_iter = self.trainer.iter
else:
self.iter += 1
if self.iter == self.patience:
+ # Early stopping condition met
self.trainer.early_stop = True
print("Early stopping occurs in iter {}, max ap is {}".format(self.best_iter, self.max_ap))
self.trainer.storage.put_scalars(timetest=12)
def after_train(self):
+ """
+ Hook to be called after training is complete to load the best model checkpoint based on AP50.
+
+ - Selects and loads the model checkpoint with the best AP50.
+ """
+ if not self.trainer.APs:
+ print("No APs were recorded during training. Skipping model selection.")
+ return
# Select the model with the best AP50
index = self.trainer.APs.index(max(self.trainer.APs)) + 1
- # Error in demo:
- # AssertionError: Checkpoint /__w/detectree2/detectree2/detectree2-data/paracou-out/train_outputs-1/model_1.pth
- # not found!
- # Therefore sleep is attempt to allow CI to pass, but it often still fails.
+ # Error handling for checkpoint loading, with a sleep to ensure file availability in CI environments
time.sleep(15)
self.trainer.checkpointer.load(self.trainer.cfg.OUTPUT_DIR + '/model_' + str(index) + '.pth')
# See https://jss367.github.io/data-augmentation-in-detectron2.html for data augmentation advice
class MyTrainer(DefaultTrainer):
- """Summary.
+ """
+ Custom Trainer class that extends the DefaultTrainer.
- Args:
- DefaultTrainer (_type_): _description_
+ This trainer adds flexibility for handling different image types (e.g., RGB and multi-band images)
+ and custom training behavior, such as early stopping and specialized data augmentation strategies.
- Returns:
- _type_: _description_
+ Args:
+ cfg (CfgNode): Configuration object containing the model and dataset configurations.
+ patience (int): Number of evaluation periods to wait for improvement before early stopping.
"""
def __init__(self, cfg, patience): # noqa: D107
self.patience = patience
- # self.resize = resize
super().__init__(cfg)
def train(self):
- """Run training.
+ """
+ Run the training loop.
- Args:
- start_iter, max_iter (int): See docs above
+ This method overrides the DefaultTrainer's train method to include early stopping and
+ custom logging of Average Precision (AP) metrics.
Returns:
- OrderedDict of results, if evaluation is enabled. Otherwise None.
+ OrderedDict: Results from evaluation, if evaluation is enabled. Otherwise, None.
"""
start_iter = self.start_iter
@@ -219,6 +368,7 @@ def train(self):
raise
finally:
self.after_train()
+ # Verify the results if testing is enabled and this is the main process
if len(self.cfg.TEST.EXPECTED_RESULTS) and comm.is_main_process():
assert hasattr(self, "_last_eval_results"), "No evaluation results obtained during training!"
verify_results(self.cfg, self._last_eval_results)
@@ -226,16 +376,61 @@ def train(self):
@classmethod
def build_evaluator(cls, cfg, dataset_name, output_folder=None):
+ """
+ Build the evaluator for the model.
+
+ Args:
+ cfg (CfgNode): Configuration object.
+ dataset_name (str): Name of the dataset to evaluate.
+ output_folder (str, optional): Directory to save evaluation results. Defaults to "eval".
+
+ Returns:
+ COCOEvaluator: An evaluator for COCO-style datasets.
+ """
if output_folder is None:
os.makedirs("eval", exist_ok=True)
output_folder = "eval"
return COCOEvaluator(dataset_name, cfg, True, output_folder)
def build_hooks(self):
+ """
+ Build the training hooks, including the custom LossEvalHook.
+
+ This method adds a custom hook for evaluating the model's loss during training, with support for
+ early stopping based on the AP50 metric.
+
+ Returns:
+ list: A list of hooks to be used during training.
+ """
hooks = super().build_hooks()
- # augmentations = [T.ResizeShortestEdge(short_edge_length=(1000, 1000),
- # max_size=1333,
- # sample_style='choice')]
+
+ # Determine the appropriate resize strategy based on the configuration
+ if self.cfg.RESIZE == "random":
+ size = None
+ # Attempt to determine the image size from the training dataset
+ for i, datas in enumerate(DatasetCatalog.get(self.cfg.DATASETS.TRAIN[0])):
+ location = datas['file_name']
+ try:
+ # Attempt to read the image with OpenCV (for RGB images)
+ img = cv2.imread(location)
+ if img is not None:
+ size = img.shape[0]
+ else:
+ # Fall back to rasterio for multi-band images
+ with rasterio.open(location) as src:
+ size = src.height # Assuming square images
+ except Exception as e:
+ # Handle any errors that occur during loading
+ print(f"Error loading image {location}: {e}")
+ continue
+ break
+ # Define augmentation based on the determined size
+ augmentations = [T.ResizeShortestEdge([size, size], size + 300)]
+ else:
+ # Use fixed size resizing as a default
+ augmentations = [T.ResizeShortestEdge([1000, 1000], 1333)]
+
+ # Insert the custom LossEvalHook before the last hook (typically the evaluation hook)
hooks.insert(
-1,
LossEvalHook(
@@ -244,53 +439,102 @@ def build_hooks(self):
build_detection_test_loader(
self.cfg,
self.cfg.DATASETS.TEST,
- DatasetMapper(self.cfg, True)
+ FlexibleDatasetMapper(self.cfg, True, augmentations=augmentations)
),
self.patience,
),
)
return hooks
+ @classmethod
+ def build_train_loader(cls, cfg):
+ """
+ Build the training data loader with support for custom augmentations and image types.
-def build_train_loader(cls, cfg):
- """Summary.
+ This method configures the data loader to apply specific augmentations depending on the image mode
+ (RGB or multi-band) and resize strategy defined in the configuration.
- Args:
- cfg (_type_): _description_
+ Args:
+ cfg (CfgNode): Configuration object.
- Returns:
- _type_: _description_
- """
- augmentations = [
- T.RandomBrightness(0.8, 1.8),
- T.RandomContrast(0.6, 1.3),
- T.RandomSaturation(0.8, 1.4),
- T.RandomRotation(angle=[90, 90], expand=False),
- T.RandomLighting(0.7),
- T.RandomFlip(prob=0.4, horizontal=True, vertical=False),
- T.RandomFlip(prob=0.4, horizontal=False, vertical=True),
- ]
+ Returns:
+ DataLoader: A data loader for the training dataset.
+ """
- if cfg.RESIZE:
- augmentations.append(T.Resize((1000, 1000)))
- elif cfg.RESIZE == "random":
- for i, datas in enumerate(DatasetCatalog.get(cfg.DATASETS.TRAIN[0])):
- location = datas['file_name']
- size = cv2.imread(location).shape[0]
- break
- print("ADD RANDOM RESIZE WITH SIZE = ", size)
- augmentations.append(T.ResizeScale(0.6, 1.4, size, size))
- return build_detection_train_loader(
- cfg,
- mapper=DatasetMapper(
+ # Define basic augmentations including rotation and flipping
+ augmentations = [
+ T.RandomRotation(angle=[90, 90], expand=False),
+ T.RandomFlip(prob=0.4, horizontal=True, vertical=False),
+ T.RandomFlip(prob=0.4, horizontal=False, vertical=True),
+ ]
+
+ # Additional augmentations for RGB images
+ if cfg.IMGMODE == "rgb":
+ augmentations.extend([
+ T.RandomBrightness(0.7, 1.5),
+ T.RandomLighting(0.7),
+ T.RandomContrast(0.6, 1.3),
+ T.RandomSaturation(0.8, 1.4)
+ ])
+
+ # Add resizing augmentations based on the resize strategy
+ if cfg.RESIZE == "fixed":
+ augmentations.append(T.ResizeShortestEdge([1000, 1000], 1333))
+ elif cfg.RESIZE == "random":
+ size = None
+ for i, datas in enumerate(DatasetCatalog.get(cfg.DATASETS.TRAIN[0])):
+ location = datas['file_name']
+ try:
+ # Try to read with cv2 (for RGB images)
+ img = cv2.imread(location)
+ if img is not None:
+ size = img.shape[0]
+ else:
+ # Fall back to rasterio for multi-band images
+ with rasterio.open(location) as src:
+ size = src.height # Assuming square images
+ except Exception as e:
+ # Handle any errors that occur during loading
+ print(f"Error loading image {location}: {e}")
+ continue
+ break
+
+ if size:
+ print("ADD RANDOM RESIZE WITH SIZE = ", size)
+ augmentations.append(T.ResizeScale(0.6, 1.4, size, size))
+ else:
+ raise ValueError("Failed to determine image size for random resize")
+ elif cfg.RESIZE == "rand_fixed":
+ augmentations.append(T.ResizeScale(0.6, 1.4, 1000, 1000))
+
+ return build_detection_train_loader(
cfg,
- is_train=True,
- augmentations=augmentations,
- ),
- )
+ mapper=FlexibleDatasetMapper(
+ cfg,
+ is_train=True,
+ augmentations=augmentations,
+ ),
+ )
+
+ @classmethod
+ def build_test_loader(cls, cfg, dataset_name):
+ """
+ Build the test data loader.
+
+ This method configures the data loader for evaluation, using the FlexibleDatasetMapper
+ to handle custom augmentations and image types.
+
+ Args:
+ cfg (CfgNode): Configuration object.
+ dataset_name (str): Name of the dataset to load for testing.
+
+ Returns:
+ DataLoader: A data loader for the test dataset.
+ """
+ return build_detection_test_loader(cfg, dataset_name, mapper=FlexibleDatasetMapper(cfg, is_train=False))
-def get_tree_dicts(directory: str, classes: List[str] = None, classes_at: str = None) -> List[Dict]:
+def get_tree_dicts(directory: str, class_mapping: Optional[Dict[str, int]] = None) -> List[Dict[str, Any]]:
"""Get the tree dictionaries.
Args:
@@ -302,36 +546,24 @@ def get_tree_dicts(directory: str, classes: List[str] = None, classes_at: str =
List of dictionaries corresponding to segmentations of trees. Each dictionary includes
bounding box around tree and points tracing a polygon around a tree.
"""
- # filepath = '/content/drive/MyDrive/forestseg/paracou_data/Panayiotis_Outputs/220303_AllSpLabelled.gpkg'
- # datagpd = gpd.read_file(filepath)
- # List_Genus = datagpd.Genus_Species.to_list()
- # Genus_Species_UniqueList = list(set(List_Genus))
-
- #
- if classes is not None:
- # list_of_classes = crowns[variable].unique().tolist()
- classes = classes
- else:
- classes = ["tree"]
- # classes = Genus_Species_UniqueList #['tree'] # genus_species list
+
dataset_dicts = []
- # for root, dirs, files in os.walk(train_location):
- # for file in files:
- # if file.endswith(".geojson"):
- # print(os.path.join(root, file))
for filename in [file for file in os.listdir(directory) if file.endswith(".geojson")]:
json_file = os.path.join(directory, filename)
with open(json_file) as f:
img_anns = json.load(f)
- # Turn off type checking for annotations until we have a better solution
- record: Dict[str, Any] = {}
- # filename = os.path.join(directory, img_anns["imagePath"])
+ record: Dict[str, Any] = {}
filename = img_anns["imagePath"]
# Make sure we have the correct height and width
- height, width = cv2.imread(filename).shape[:2]
+ # If image path ends in .png use cv2 to get height and width else if image path ends in .tif use rasterio
+ if filename.endswith(".png"):
+ height, width = cv2.imread(filename).shape[:2]
+ elif filename.endswith(".tif"):
+ with rasterio.open(filename) as src:
+ height, width = src.shape
record["file_name"] = filename
record["height"] = height
@@ -343,65 +575,78 @@ def get_tree_dicts(directory: str, classes: List[str] = None, classes_at: str =
objs = []
for features in img_anns["features"]:
anno = features["geometry"]
- # pdb.set_trace()
- # GenusSpecies = features['properties']['Genus_Species']
px = [a[0] for a in anno["coordinates"][0]]
py = [np.array(height) - a[1] for a in anno["coordinates"][0]]
- # print("### HERE IS PY ###", py)
poly = [(x, y) for x, y in zip(px, py)]
poly = [p for x in poly for p in x]
- # print("#### HERE ARE SOME POLYS #####", poly)
- if classes != ["tree"]:
- obj = {
- "bbox": [np.min(px), np.min(py), np.max(px), np.max(py)],
- "bbox_mode": BoxMode.XYXY_ABS,
- "segmentation": [poly],
- "category_id": classes.index(features["properties"][classes_at]), # id
- # "category_id": 0, #id
- "iscrowd": 0,
- }
+
+ # If class mapping is provided, use it; otherwise, default to "tree"
+ if class_mapping:
+ category_id = class_mapping[features["properties"]["status"]]
else:
- obj = {
- "bbox": [np.min(px), np.min(py), np.max(px), np.max(py)],
- "bbox_mode": BoxMode.XYXY_ABS,
- "segmentation": [poly],
- "category_id": 0, # id
- "iscrowd": 0,
- }
- # pdb.set_trace()
+ category_id = 0 # Default to "tree" if no class mapping is provided
+
+ obj = {
+ "bbox": [np.min(px), np.min(py), np.max(px), np.max(py)],
+ "bbox_mode": BoxMode.XYXY_ABS,
+ "segmentation": [poly],
+ "category_id": category_id,
+ "iscrowd": 0,
+ }
+
objs.append(obj)
- # print("#### HERE IS OBJS #####", objs)
- record["annotations"] = objs
+
+ record["annotations"] = objs if objs else []
dataset_dicts.append(record)
+
return dataset_dicts
def combine_dicts(root_dir: str,
val_dir: int,
mode: str = "train",
- classes: List[str] = None,
- classes_at: str = None) -> List[Dict]:
- """Join tree dicts from different directories.
+ class_mapping: Optional[Dict[str, int]] = None) -> List[Dict[str, Any]]:
+ """
+ Combine dictionaries from different directories based on the specified mode.
+
+ This function aggregates tree dictionaries from multiple directories within a root directory.
+ Depending on the mode, it either combines dictionaries from all directories,
+ all except a specified validation directory, or only from the validation directory.
Args:
- root_dir:
- val_dir:
+ root_dir (str): The root directory containing subdirectories with tree dictionaries.
+ val_dir (int): The index (1-based) of the validation directory to exclude or use depending on the mode.
+ mode (str, optional): The mode of operation. Can be "train", "val", or "full".
+ "train" excludes the validation directory,
+ "val" includes only the validation directory,
+ and "full" includes all directories. Defaults to "train".
+ class_mapping: A dictionary mapping class labels to category indices (optional).
Returns:
- Concatenated array of dictionaries over all directories
+ List of combined dictionaries from the specified directories.
"""
- train_dirs = [os.path.join(root_dir, dir) for dir in os.listdir(root_dir)]
+ # Get the list of directories within the root directory
+ train_dirs = [
+ os.path.join(root_dir, dir)
+ for dir in os.listdir(root_dir)
+ if os.path.isdir(os.path.join(root_dir, dir))
+ ]
+ # Handle the different modes for combining dictionaries
if mode == "train":
+ # Exclude the validation directory from the list of directories
del train_dirs[(val_dir - 1)]
tree_dicts = []
for d in train_dirs:
- tree_dicts += get_tree_dicts(d, classes=classes, classes_at=classes_at)
+ # Combine dictionaries from all directories except the validation directory
+ tree_dicts += get_tree_dicts(d, class_mapping=class_mapping)
elif mode == "val":
- tree_dicts = get_tree_dicts(train_dirs[(val_dir - 1)], classes=classes, classes_at=classes_at)
+ # Use only the validation directory
+ tree_dicts = get_tree_dicts(train_dirs[(val_dir - 1)], class_mapping=class_mapping)
elif mode == "full":
+ # Combine dictionaries from all directories, including the validation directory
tree_dicts = []
for d in train_dirs:
- tree_dicts += get_tree_dicts(d, classes=classes, classes_at=classes_at)
+ tree_dicts += get_tree_dicts(d, class_mapping=class_mapping)
return tree_dicts
@@ -427,41 +672,52 @@ def get_filenames(directory: str):
def register_train_data(train_location,
name: str = "tree",
val_fold=None,
- classes=None,
- classes_at=None):
+ class_mapping_file=None):
"""Register data for training and (optionally) validation.
Args:
- train_location: directory containing training folds
- name: string to name data
- val_fold: fold assigned for validation and tuning. If not given,
- will take place on all folds.
+ train_location: Directory containing training folds.
+ name: Name to register the dataset.
+ val_fold: Validation fold index (optional).
+ class_mapping_file: Path to the class mapping file (json or pickle).
"""
+ # Load the class mapping from file if provided
+ class_mapping = None
+ if class_mapping_file:
+ class_mapping = load_class_mapping(class_mapping_file)
+ thing_classes = list(class_mapping.keys()) # Convert dictionary to list of class names
+ print(f"Class mapping loaded: {class_mapping}") # Debugging step
+ else:
+ thing_classes = ["tree"]
+
if val_fold is not None:
for d in ["train", "val"]:
- DatasetCatalog.register(name + "_" + d, lambda d=d: combine_dicts(train_location,
- val_fold, d,
- classes=classes, classes_at=classes_at))
- if classes is None:
- MetadataCatalog.get(name + "_" + d).set(thing_classes=["tree"])
- else:
- MetadataCatalog.get(name + "_" + d).set(thing_classes=classes)
+ DatasetCatalog.register(
+ name + "_" + d,
+ lambda d=d: combine_dicts(train_location, val_fold, d, class_mapping=class_mapping)
+ )
+ MetadataCatalog.get(name + "_" + d).set(thing_classes=thing_classes)
else:
- DatasetCatalog.register(name + "_" + "full", lambda d=d: combine_dicts(train_location,
- 0, "full",
- classes=classes, classes_at=classes_at))
- if classes is None:
- MetadataCatalog.get(name + "_" + "full").set(thing_classes=["tree"])
- else:
- MetadataCatalog.get(name + "_" + "full").set(thing_classes=classes)
+ DatasetCatalog.register(
+ name + "_" + "full",
+ lambda d=d: combine_dicts(train_location, 0, "full", class_mapping=class_mapping)
+ )
+ MetadataCatalog.get(name + "_" + "full").set(thing_classes=thing_classes)
+
+def get_classes(out_dir):
+ """Function that will read the classes that are recorded during tiling.
-def read_data(out_dir):
- """Function that will read the classes that are recorded during tiling."""
+ Args:
+ out_dir: directory where classes.txt is located
+
+ Returns:
+ list of classes
+ """
list = []
- out_tif = out_dir + 'classes.txt'
+ classes_txt = out_dir + 'classes.txt'
# open file and read the content in a list
- with open(out_tif, 'r') as fp:
+ with open(classes_txt, 'r') as fp:
for line in fp:
# remove linebreak from a current name
# linebreak is the last character of each line
@@ -483,14 +739,32 @@ def remove_registered_data(name="tree"):
def register_test_data(test_location, name="tree"):
- """Register data for testing."""
+ """Register data for testing.
+
+ Args:
+ test_location: directory containing test data
+ name: string to name data
+ """
d = "test"
- DatasetCatalog.register(name + "_" + d, lambda d=d: get_tree_dicts(test_location))
- MetadataCatalog.get(name + "_" + d).set(thing_classes=["tree"])
+
+ class_mapping = None
+ if class_mapping_file:
+ class_mapping = load_class_mapping(class_mapping_file)
+ thing_classes = list(class_mapping.keys()) # Convert dictionary to list of class names
+ print(f"Class mapping loaded: {class_mapping}") # Debugging step
+ else:
+ thing_classes = ["tree"]
+
+ DatasetCatalog.register(name + "_" + d, lambda d=d: get_tree_dicts(test_location, class_mapping))
+ MetadataCatalog.get(name + "_" + d).set(thing_classes=thing_classes)
def load_json_arr(json_path):
- """Load json array."""
+ """Load json array.
+
+ Args:
+ json_path: path to json file
+ """
lines = []
with open(json_path, "r") as f:
for line in f:
@@ -513,10 +787,12 @@ def setup_cfg(
base_lr=0.0003389,
weight_decay=0.001,
max_iter=1000,
- num_classes=1,
eval_period=100,
out_dir="./train_outputs",
- resize=True,
+ resize="fixed", # "fixed" or "random" or "rand_fixed"
+ imgmode="rgb",
+ num_bands=3,
+ class_mapping_file=None,
):
"""Set up config object # noqa: D417.
@@ -538,7 +814,23 @@ def setup_cfg(
num_classes: number of classes
eval_period: number of iterations between evaluations
out_dir: directory to save outputs
+ resize: resize strategy for images
+ imgmode: image mode (rgb or multispectral)
+ num_bands: number of bands in the image
+ class_mapping_file: path to class mapping file
"""
+
+ # Load the class mapping if provided
+ if class_mapping_file:
+ class_mapping = load_class_mapping(class_mapping_file)
+ num_classes = len(class_mapping) # Set the number of classes based on the mapping
+ else:
+ num_classes = 1 # Default to 1 class if no mapping is provided
+
+ # Validate the resize parameter
+ if resize not in {"fixed", "random", "rand_fixed"}:
+ raise ValueError(f"Invalid resize option '{resize}'. Must be 'fixed', 'random', or 'rand_fixed'.")
+
cfg = get_cfg()
cfg.merge_from_file(model_zoo.get_config_file(base_model))
cfg.DATASETS.TRAIN = trains
@@ -566,6 +858,20 @@ def setup_cfg(
cfg.TEST.EVAL_PERIOD = eval_period
cfg.RESIZE = resize
cfg.INPUT.MIN_SIZE_TRAIN = 1000
+ cfg.IMGMODE = imgmode # "rgb" or "ms" (multispectral)
+ if num_bands > 3:
+ # Adjust PIXEL_MEAN and PIXEL_STD for the number of bands
+ default_pixel_mean = cfg.MODEL.PIXEL_MEAN
+ default_pixel_std = cfg.MODEL.PIXEL_STD
+ # Extend or truncate the PIXEL_MEAN and PIXEL_STD based on num_bands
+ cfg.MODEL.PIXEL_MEAN = (
+ default_pixel_mean * (num_bands // len(default_pixel_mean))
+ + default_pixel_mean[:num_bands % len(default_pixel_mean)]
+ )
+ cfg.MODEL.PIXEL_STD = (
+ default_pixel_std * (num_bands // len(default_pixel_std))
+ + default_pixel_std[:num_bands % len(default_pixel_std)]
+ )
return cfg
@@ -576,7 +882,17 @@ def predictions_on_data(directory=None,
scale=1,
geos_exist=True,
num_predictions=0):
- """Prediction produced from a test folder and outputted to predictions folder."""
+ """Prediction produced from a test folder and outputted to predictions folder.
+
+ Args:
+ directory: directory containing test data
+ predictor: predictor object
+ trees_metadata: metadata for trees
+ save: boolean to save predictions
+ scale: scale of image
+ geos_exist: boolean to determine if geojson files exist
+ num_predictions: number of predictions to make
+ """
test_location = directory + "/test"
pred_dir = test_location + "/predictions"
@@ -623,46 +939,94 @@ def predictions_on_data(directory=None,
json.dump(evaluations, dest)
+def modify_conv1_weights(model, num_input_channels):
+ """
+ Modify the weights of the first convolutional layer (conv1) to accommodate a different number of input channels.
+
+ This function adjusts the weights of the `conv1` layer in the model's backbone to support a custom number
+ of input channels. It creates a new weight tensor with the desired number of input channels,
+ and initializes it by repeating the weights of the original channels.
+
+ Args:
+ model (torch.nn.Module): The model containing the convolutional layer to modify.
+ num_input_channels (int): The number of input channels for the new conv1 layer.
+
+ """
+ with torch.no_grad():
+ # Retrieve the original weights of the conv1 layer
+ old_weights = model.backbone.bottom_up.stem.conv1.weight
+
+ # Create a new weight tensor with the desired number of input channels
+ # The shape is (out_channels, in_channels, height, width)
+ new_weights = torch.zeros((old_weights.size(0), num_input_channels, *old_weights.shape[2:]))
+
+ # Initialize the new weights by repeating the original weights across the new channels
+ # This example repeats the first 3 channels if num_input_channels > 3
+ for i in range(num_input_channels):
+ new_weights[:, i, :, :] = old_weights[:, i % 3, :, :]
+
+ # Create a new conv1 layer with the updated number of input channels
+ model.backbone.bottom_up.stem.conv1 = nn.Conv2d(
+ num_input_channels, old_weights.size(0), kernel_size=7, stride=2, padding=3, bias=False
+ )
+
+ # Copy the modified weights into the new conv1 layer
+ model.backbone.bottom_up.stem.conv1.weight.copy_(new_weights)
+
+
+def get_latest_model_path(output_dir: str) -> str:
+ """
+ Find the model file with the highest index in the specified output directory.
+
+ Args:
+ output_dir (str): The directory where the model files are stored.
+
+ Returns:
+ str: The path to the model file with the highest index.
+ """
+ # Regular expression to match model files with the pattern "model_X.pth"
+ model_pattern = re.compile(r"model_(\d+)\.pth")
+
+ # List all files in the output directory
+ files = os.listdir(output_dir)
+
+ # Find all files that match the pattern and extract their indices
+ model_files = []
+ for f in files:
+ match = model_pattern.search(f)
+ if match:
+ model_files.append((f, int(match.group(1))))
+
+ if not model_files:
+ raise FileNotFoundError(f"No model files found in the directory {output_dir}")
+
+ # Sort the files by index in descending order and select the highest one
+ latest_model_file = max(model_files, key=lambda x: x[1])[0]
+
+ # Return the full path to the latest model file
+ return os.path.join(output_dir, latest_model_file)
+
+
if __name__ == "__main__":
- train_location = "/content/drive/Shareddrives/detectree2/data/Paracou/tiles/train/"
- register_train_data(train_location, "Paracou", 1) # folder, name, validation fold
-
- name = "Paracou2019"
- train_location = "/content/drive/Shareddrives/detectree2/data/Paracou/tiles2019/train/"
- dataset_dicts = combine_dicts(train_location, 1)
- trees_metadata = MetadataCatalog.get(name + "_train")
- # dataset_dicts = get_tree_dicts("./")
- for d in dataset_dicts:
- img = cv2.imread(d["file_name"])
- visualizer = Visualizer(img[:, :, ::-1], metadata=trees_metadata, scale=0.5)
- out = visualizer.draw_dataset_dict(d)
- image = cv2.cvtColor(out.get_image()[:, :, ::-1], cv2.COLOR_BGR2RGB)
- # display(Image.fromarray(image))
- # Set the base (pre-trained) model from the detectron2 model_zoo
- model = "COCO-InstanceSegmentation/mask_rcnn_R_101_FPN_3x.yaml"
- # Set the names of the registered train and test sets
- # pretrained model?
- # trained_model = "/content/drive/Shareddrives/detectree2/models/220629_ParacouSepilokDanum_JB.pth"
- trains = (
- "Paracou_train",
- "Paracou2019_train",
- "ParacouUAV_train",
- "Danum_train",
- "SepilokEast_train",
- "SepilokWest_train",
- )
- tests = (
- "Paracou_val",
- "Paracou2019_val",
- "ParacouUAV_val",
- "Danum_val",
- "SepilokEast_val",
- "SepilokWest_val",
+ # Define paths to training data and optional class mapping file
+ train_location = "/path/to/your/train/location"
+ class_mapping_file = "/path/to/your/class_to_idx.json" # Optional, can be None
+
+ # Register the training and validation datasets using the class mapping
+ # If class_mapping_file is not provided, defaults to "tree"
+ register_train_data(train_location, "MyDataset", val_fold=1, class_mapping_file=class_mapping_file)
+
+ # Set up model configuration, using the class mapping to determine the number of classes
+ cfg = setup_cfg(
+ base_model="COCO-InstanceSegmentation/mask_rcnn_R_101_FPN_3x.yaml",
+ trains=("MyDataset_train", ),
+ tests=("MyDataset_val", ),
+ max_iter=3000,
+ out_dir="/path/to/output",
+ class_mapping_file=class_mapping_file # Optional
)
- out_dir = "/content/drive/Shareddrives/detectree2/220703_train_outputs"
- # update_model arg can be used to load in trained model
- cfg = setup_cfg(model, trains, tests, eval_period=100, max_iter=3000, out_dir=out_dir)
+ # Train the model
trainer = MyTrainer(cfg, patience=4)
trainer.resume_or_load(resume=False)
trainer.train()
diff --git a/detectree2/preprocessing/tiling.py b/detectree2/preprocessing/tiling.py
index 2ceeb8c4..2a9c1d7d 100644
--- a/detectree2/preprocessing/tiling.py
+++ b/detectree2/preprocessing/tiling.py
@@ -4,11 +4,14 @@
of models and making landscape predictions.
"""
+import concurrent.futures
import json
+import logging
import os
+import pickle
import random
import shutil
-import warnings
+import warnings # noqa: F401
from math import ceil
from pathlib import Path
@@ -17,11 +20,17 @@
import numpy as np
import rasterio
from fiona.crs import from_epsg # noqa: F401
-from rasterio.crs import CRS
-from rasterio.io import DatasetReader
+# from rasterio.crs import CRS
+from rasterio.errors import RasterioIOError
+# from rasterio.io import DatasetReader
from rasterio.mask import mask
+# from rasterio.windows import from_bounds
from shapely.geometry import box
+# Configure logging
+logging.basicConfig(level=logging.INFO)
+logger = logging.getLogger(__name__)
+
# class img_data(DatasetReader):
# """
# Class for image data to be processed for tiling
@@ -47,87 +56,93 @@ def get_features(gdf: gpd.GeoDataFrame):
return [json.loads(gdf.to_json())["features"][0]["geometry"]]
-def tile_data(
- data: DatasetReader,
- out_dir: str,
- buffer: int = 30,
- tile_width: int = 200,
- tile_height: int = 200,
- dtype_bool: bool = False,
-) -> None:
- """Tiles up orthomosaic for making predictions on.
+def load_class_mapping(file_path: str):
+ """Function to load class-to-index mapping from a file.
- Tiles up full othomosaic into managable chunks to make predictions on. Use tile_data_train to generate tiled
- training data. A bug exists on some input raster types whereby outputed tiles are completely black - the dtype_bool
- argument should be switched if this is the case.
+ Args:
+ file_path: Path to the file (json or pickle)
+
+ Returns:
+ class_to_idx: Loaded class-to-index mapping
+ """
+ file_ext = Path(file_path).suffix
+
+ if file_ext == '.json':
+ with open(file_path, 'r') as f:
+ class_to_idx = json.load(f)
+ elif file_ext == '.pkl':
+ with open(file_path, 'rb') as f:
+ class_to_idx = pickle.load(f)
+ else:
+ raise ValueError("Unsupported file format. Use '.json' or '.pkl'.")
+
+ return class_to_idx
+
+
+def process_tile(
+ img_path: str,
+ out_dir: str,
+ buffer: int,
+ tile_width: int,
+ tile_height: int,
+ dtype_bool: bool,
+ minx,
+ miny,
+ crs,
+ tilename,
+ crowns: gpd.GeoDataFrame = None,
+ threshold: float = 0,
+ nan_threshold: float = 0,
+):
+ """Process a single tile for making predictions.
Args:
- data: Orthomosaic as a rasterio object in a UTM type projection
+ img_path: Path to the orthomosaic
+ out_dir: Output directory
buffer: Overlapping buffer of tiles in meters (UTM)
tile_width: Tile width in meters
tile_height: Tile height in meters
dtype_bool: Flag to edit dtype to prevent black tiles
+ minx: Minimum x coordinate of tile
+ miny: Minimum y coordinate of tile
+ crs: Coordinate reference system
+ tilename: Name of the tile
Returns:
None
"""
- out_path = Path(out_dir)
- os.makedirs(out_path, exist_ok=True)
- crs = CRS.from_string(data.crs.wkt)
- crs = crs.to_epsg()
- tilename = Path(data.name).stem
-
- total_tiles = int(
- ((data.bounds[2] - data.bounds[0]) / tile_width) * ((data.bounds[3] - data.bounds[1]) / tile_height)
- )
-
- tile_count = 0
- print(f"Tiling to {total_tiles} total tiles")
-
- for minx in np.arange(data.bounds[0], data.bounds[2] - tile_width,
- tile_width, int):
- for miny in np.arange(data.bounds[1], data.bounds[3] - tile_height,
- tile_height, int):
-
- tile_count += 1
- # Naming conventions
+ try:
+ with rasterio.open(img_path) as data:
+ out_path = Path(out_dir)
out_path_root = out_path / f"{tilename}_{minx}_{miny}_{tile_width}_{buffer}_{crs}"
- # new tiling bbox including the buffer
- bbox = box(
- minx - buffer,
- miny - buffer,
- minx + tile_width + buffer,
- miny + tile_height + buffer,
- )
- # define the bounding box of the tile, excluding the buffer
- # (hence selecting just the central part of the tile)
- # bbox_central = box(minx, miny, minx + tile_width, miny + tile_height)
- # turn the bounding boxes into geopandas DataFrames
- geo = gpd.GeoDataFrame({"geometry": bbox}, index=[0], crs=data.crs)
- # geo_central = gpd.GeoDataFrame(
- # {"geometry": bbox_central}, index=[0], crs=from_epsg(4326)
- # ) # 3182
- # overlapping_crowns = sjoin(crowns, geo_central, how="inner")
+ minx_buffered = minx - buffer
+ miny_buffered = miny - buffer
+ maxx_buffered = minx + tile_width + buffer
+ maxy_buffered = miny + tile_height + buffer
- # here we are cropping the tiff to the bounding box of the tile we want
+ bbox = box(minx_buffered, miny_buffered, maxx_buffered, maxy_buffered)
+ geo = gpd.GeoDataFrame({"geometry": bbox}, index=[0], crs=data.crs)
coords = get_features(geo)
- # print("Coords:", coords)
- # define the tile as a mask of the whole tiff with just the bounding box
+ overlapping_crowns = None
+ if crowns is not None:
+ overlapping_crowns = gpd.clip(crowns, geo)
+ if overlapping_crowns.empty or (overlapping_crowns.dissolve().area[0] / geo.area[0]) < threshold:
+ return None
+
out_img, out_transform = mask(data, shapes=coords, crop=True)
- # Discard scenes with many out-of-range pixels
- out_sumbands = np.sum(out_img, 0)
+ out_sumbands = np.sum(out_img, axis=0)
zero_mask = np.where(out_sumbands == 0, 1, 0)
nan_mask = np.where(out_sumbands == 765, 1, 0)
sumzero = zero_mask.sum()
sumnan = nan_mask.sum()
totalpix = out_img.shape[1] * out_img.shape[2]
- if sumzero > 0.25 * totalpix:
- continue
- elif sumnan > 0.25 * totalpix:
- continue
+
+ # If the tile is mostly empty or mostly nan, don't save it
+ if sumzero > nan_threshold * totalpix or sumnan > nan_threshold * totalpix:
+ return None
out_meta = data.meta.copy()
out_meta.update({
@@ -137,145 +152,103 @@ def tile_data(
"transform": out_transform,
"nodata": None,
})
- # dtype needs to be unchanged for some data and set to uint8 for others
if dtype_bool:
out_meta.update({"dtype": "uint8"})
- # print("Out Meta:",out_meta)
- # Saving the tile as a new tiff, named by the origin of the tile.
- # If tile appears blank in folder can show the image here and may
- # need to fix RGB data or the dtype
- # show(out_img)
- out_tif = out_path_root.with_suffix(out_path_root.suffix + ".tif")
+ out_tif = out_path_root.with_suffix(".tif")
with rasterio.open(out_tif, "w", **out_meta) as dest:
dest.write(out_img)
- # read in the tile we have just saved
- clipped = rasterio.open(out_tif)
- # read it as an array
- # show(clipped)
- arr = clipped.read()
+ with rasterio.open(out_tif) as clipped:
+ arr = clipped.read()
+ r, g, b = arr[0], arr[1], arr[2]
+ rgb = np.dstack((b, g, r)) # Reorder for cv2 (BGRA)
- # each band of the tiled tiff is a colour!
- r = arr[0]
- g = arr[1]
- b = arr[2]
+ # Rescale to 0-255 if necessary
+ if np.max(g) > 255:
+ rgb_rescaled = 255 * rgb / 65535
+ else:
+ rgb_rescaled = rgb
- # stack up the bands in an order appropriate for saving with cv2,
- # then rescale to the correct 0-255 range for cv2
+ cv2.imwrite(str(out_path_root.with_suffix(".png").resolve()), rgb_rescaled)
- rgb = np.dstack((b, g, r)) # BGR for cv2
+ if overlapping_crowns is not None:
+ return data, out_path_root, overlapping_crowns, minx, miny, buffer
- if np.max(g) > 255:
- rgb_rescaled = 255 * rgb / 65535
- else:
- rgb_rescaled = rgb # scale to image
- # print("rgb rescaled", rgb_rescaled)
-
- # save this as jpg or png...we are going for png...again, named with the origin of the specific tile
- # here as a naughty method
- cv2.imwrite(
- str(out_path_root.with_suffix(out_path_root.suffix + ".png").resolve()),
- rgb_rescaled,
- )
- if tile_count % 50 == 0:
- print(f"Processed {tile_count} tiles of {total_tiles} tiles")
+ return data, out_path_root, None, minx, miny, buffer
- print("Tiling complete")
+ except RasterioIOError as e:
+ logger.error(f"RasterioIOError while applying mask {coords}: {e}")
+ return None
+ except Exception as e:
+ logger.error(f"Error processing tile {tilename} at ({minx}, {miny}): {e}")
+ return None
-def tile_data_train( # noqa: C901
- data: DatasetReader,
+def process_tile_ms(
+ img_path: str,
out_dir: str,
- buffer: int = 30,
- tile_width: int = 200,
- tile_height: int = 200,
+ buffer: int,
+ tile_width: int,
+ tile_height: int,
+ dtype_bool: bool,
+ minx,
+ miny,
+ crs,
+ tilename,
crowns: gpd.GeoDataFrame = None,
threshold: float = 0,
- nan_threshold: float = 0.1,
- dtype_bool: bool = False,
-) -> None:
- """Tiles up orthomosaic and corresponding crowns into training tiles.
-
- A threshold can be used to ensure a good coverage of crowns across a tile. Tiles that do not have sufficient
- coverage are rejected.
+ nan_threshold: float = 0,
+):
+ """Process a single tile for making predictions.
Args:
- data: Orthomosaic as a rasterio object in a UTM type projection
+ img_path: Path to the orthomosaic
+ out_dir: Output directory
buffer: Overlapping buffer of tiles in meters (UTM)
tile_width: Tile width in meters
tile_height: Tile height in meters
- crowns: Crown polygons as a geopandas dataframe
- threshold: Min proportion of the tile covered by crowns to be accepted {0,1}
- nan_theshold: Max proportion of tile covered by nans
dtype_bool: Flag to edit dtype to prevent black tiles
+ minx: Minimum x coordinate of tile
+ miny: Minimum y coordinate of tile
+ crs: Coordinate reference system
+ tilename: Name of the tile
Returns:
None
-
"""
-
- # TODO: Clip data to crowns straight away to speed things up
- # TODO: Tighten up epsg handling
- out_path = Path(out_dir)
- os.makedirs(out_path, exist_ok=True)
- tilename = Path(data.name).stem
- crs = CRS.from_string(data.crs.wkt)
- crs = crs.to_epsg()
- # out_img, out_transform = mask(data, shapes=crowns.buffer(buffer), crop=True)
- # Should start from data.bounds[0] + buffer, data.bounds[1] + buffer to avoid later complications
- for minx in np.arange(ceil(data.bounds[0]) + buffer, data.bounds[2] - tile_width - buffer, tile_width, int):
- for miny in np.arange(ceil(data.bounds[1]) + buffer, data.bounds[3] - tile_height - buffer, tile_height, int):
-
+ try:
+ with rasterio.open(img_path) as data:
+ out_path = Path(out_dir)
out_path_root = out_path / f"{tilename}_{minx}_{miny}_{tile_width}_{buffer}_{crs}"
- # Calculate the buffered tile dimensions
- # tile_width_buffered = tile_width + 2 * buffer
- # tile_height_buffered = tile_height + 2 * buffer
-
- # Calculate the bounding box coordinates with buffer
minx_buffered = minx - buffer
miny_buffered = miny - buffer
maxx_buffered = minx + tile_width + buffer
maxy_buffered = miny + tile_height + buffer
- # Create the affine transformation matrix for the tile
- # transform = from_bounds(minx_buffered, miny_buffered, maxx_buffered,
- # maxy_buffered, tile_width_buffered, tile_height_buffered)
-
bbox = box(minx_buffered, miny_buffered, maxx_buffered, maxy_buffered)
- geo = gpd.GeoDataFrame({"geometry": bbox}, index=[0], crs=data.crs)
- coords = get_features(geo)
+ geo = gpd.GeoDataFrame({"geometry": [bbox]}, index=[0], crs=data.crs)
+ coords = [geo.geometry[0].__geo_interface__]
- # Skip if insufficient coverage of crowns - good to have early on to save on unnecessary processing
- with warnings.catch_warnings():
- warnings.simplefilter("ignore")
- # Warning:
- # _crs_mismatch_warn
+ overlapping_crowns = None
+ if crowns is not None:
overlapping_crowns = gpd.clip(crowns, geo)
+ if overlapping_crowns.empty or (overlapping_crowns.dissolve().area[0] / geo.area[0]) < threshold:
+ return None
- # Ignore tiles with no crowns
- if overlapping_crowns.empty:
- continue
-
- # Discard tiles that do not have a sufficient coverage of training crowns
- if (overlapping_crowns.dissolve().area[0] / geo.area[0]) < threshold:
- continue
-
- # define the tile as a mask of the whole tiff with just the bounding box
out_img, out_transform = mask(data, shapes=coords, crop=True)
- # Discard scenes with many out-of-range pixels
- out_sumbands = np.sum(out_img, 0)
+ out_sumbands = np.sum(out_img, axis=0)
zero_mask = np.where(out_sumbands == 0, 1, 0)
- nan_mask = np.where(out_sumbands == 765, 1, 0)
+ nan_mask = np.isnan(out_sumbands)
sumzero = zero_mask.sum()
sumnan = nan_mask.sum()
totalpix = out_img.shape[1] * out_img.shape[2]
- if sumzero > nan_threshold * totalpix: # reject tiles with many 0 cells
- continue
- elif sumnan > nan_threshold * totalpix: # reject tiles with many NaN cells
- continue
+
+ # If the tile is mostly empty or mostly nan, don't save it
+ if sumzero > nan_threshold * totalpix or sumnan > nan_threshold * totalpix:
+ return None
out_meta = data.meta.copy()
out_meta.update({
@@ -285,104 +258,171 @@ def tile_data_train( # noqa: C901
"transform": out_transform,
"nodata": None,
})
- # dtype needs to be unchanged for some data and set to uint8 for others to deal with black tiles
if dtype_bool:
out_meta.update({"dtype": "uint8"})
- # Saving the tile as a new tiff, named by the origin of the tile. If tile appears blank in folder can show
- # the image here and may need to fix RGB data or the dtype
- out_tif = out_path_root.with_suffix(out_path_root.suffix + ".tif")
+ out_tif = out_path_root.with_suffix(".tif")
with rasterio.open(out_tif, "w", **out_meta) as dest:
dest.write(out_img)
- # read in the tile we have just saved
- clipped = rasterio.open(out_tif)
+ if overlapping_crowns is not None:
+ return data, out_path_root, overlapping_crowns, minx, miny, buffer
- # read it as an array
- arr = clipped.read()
+ return data, out_path_root, None, minx, miny, buffer
- # each band of the tiled tiff is a colour!
- r = arr[0]
- g = arr[1]
- b = arr[2]
+ except RasterioIOError as e:
+ logger.error(f"RasterioIOError while applying mask {coords}: {e}")
+ return None
+ except Exception as e:
+ logger.error(f"Error processing tile {tilename} at ({minx}, {miny}): {e}")
+ return None
- # stack up the bands in an order appropriate for saving with cv2, then rescale to the correct 0-255 range
- # for cv2. BGR ordering is correct for cv2 (and detectron2)
- rgb = np.dstack((b, g, r))
- # Some rasters need to have values rescaled to 0-255
- # TODO: more robust check
- if np.max(g) > 255:
- rgb_rescaled = 255 * rgb / 65535
- else:
- # scale to image
- rgb_rescaled = rgb
-
- # save this as png, named with the origin of the specific tile
- # potentially bad practice
- cv2.imwrite(
- str(out_path_root.with_suffix(out_path_root.suffix + ".png").resolve()),
- rgb_rescaled,
- )
+def process_tile_train(
+ img_path: str,
+ out_dir: str,
+ buffer: int,
+ tile_width: int,
+ tile_height: int,
+ dtype_bool: bool,
+ minx,
+ miny,
+ crs,
+ tilename,
+ crowns: gpd.GeoDataFrame,
+ threshold,
+ nan_threshold,
+ mode: str = "rgb",
+ class_column: str = None, # Allow user to specify class column
+) -> None:
+ """Process a single tile for training data.
+
+ Args:
+ img_path: Path to the orthomosaic
+ out_dir: Output directory
+ buffer: Overlapping buffer of tiles in meters (UTM)
+ tile_width: Tile width in meters
+ tile_height: Tile height in meters
+ dtype_bool: Flag to edit dtype to prevent black tiles
+ minx: Minimum x coordinate of tile
+ miny: Minimum y coordinate of tile
+ crs: Coordinate reference system
+ tilename: Name of the tile
+ crowns: Crown polygons as a geopandas dataframe
+ threshold: Min proportion of the tile covered by crowns to be accepted {0,1}
+ nan_theshold: Max proportion of tile covered by nans
+
+ Returns:
+ None
+ """
+ if mode == "rgb":
+ result = process_tile(img_path, out_dir, buffer, tile_width, tile_height, dtype_bool, minx, miny, crs,
+ tilename, crowns, threshold, nan_threshold)
+ elif mode == "ms":
+ result = process_tile_ms(img_path, out_dir, buffer, tile_width, tile_height, dtype_bool, minx, miny, crs,
+ tilename, crowns, threshold, nan_threshold)
+
+ if result is None:
+ # logger.warning(f"Skipping tile at ({minx}, {miny}) due to insufficient data.")
+ return
+
+ data, out_path_root, overlapping_crowns, minx, miny, buffer = result
+
+ overlapping_crowns = overlapping_crowns.explode(index_parts=True)
+ moved = overlapping_crowns.translate(-minx + buffer, -miny + buffer)
+ scalingx = 1 / (data.transform[0])
+ scalingy = -1 / (data.transform[4])
+ moved_scaled = moved.scale(scalingx, scalingy, origin=(0, 0))
+
+ if mode == "rgb":
+ impath = {"imagePath": out_path_root.with_suffix(".png").as_posix()}
+ elif mode == "ms":
+ impath = {"imagePath": out_path_root.with_suffix(".tif").as_posix()}
+
+ try:
+ filename = out_path_root.with_suffix(".geojson")
+ moved_scaled = overlapping_crowns.set_geometry(moved_scaled)
+
+ if class_column is not None:
+ # Ensure we map the selected column to the 'status' field
+ moved_scaled['status'] = moved_scaled[class_column]
+ # Keep only 'status' and geometry
+ moved_scaled = moved_scaled[['geometry', 'status']]
+ else:
+ # Keep only geometry to reduce file size
+ moved_scaled = moved_scaled[['geometry']]
+
+ # Save the result as GeoJSON
+ moved_scaled.to_file(driver="GeoJSON", filename=filename)
+
+ # Add image path info to the GeoJSON file
+ with open(filename, "r") as f:
+ shp = json.load(f)
+ shp.update(impath)
+ with open(filename, "w") as f:
+ json.dump(shp, f)
+ except ValueError:
+ logger.warning("Cannot write empty DataFrame to file.")
+ return
+
+
+# Define a top-level helper function
+def process_tile_train_helper(args):
+ return process_tile_train(*args)
+
+
+def tile_data(
+ img_path: str,
+ out_dir: str,
+ buffer: int = 30,
+ tile_width: int = 200,
+ tile_height: int = 200,
+ crowns: gpd.GeoDataFrame = None,
+ threshold: float = 0,
+ nan_threshold: float = 0.1,
+ dtype_bool: bool = False,
+ mode: str = "rgb",
+ class_column: str = None, # Allow class column to be passed here
+) -> None:
+ """Tiles up orthomosaic and corresponding crowns (if supplied) into training/prediction tiles.
- # select the crowns that intersect the non-buffered central
- # section of the tile using the inner join
- # TODO: A better solution would be to clip crowns to tile extent
- # overlapping_crowns = sjoin(crowns, geo_central, how="inner")
- # Maybe left join to keep information of crowns?
-
- overlapping_crowns = overlapping_crowns.explode(index_parts=True)
-
- # Translate to 0,0 to overlay on png
- moved = overlapping_crowns.translate(-minx + buffer, -miny + buffer)
-
- # scale to deal with the resolution
- scalingx = 1 / (data.transform[0])
- scalingy = -1 / (data.transform[4])
- moved_scaled = moved.scale(scalingx, scalingy, origin=(0, 0))
-
- impath = {"imagePath": out_path_root.with_suffix(out_path_root.suffix + ".png").as_posix()}
-
- # Save as a geojson, a format compatible with detectron2, again named by the origin of the tile.
- # If the box selected from the image is outside of the mapped region due to the image being on a slant
- # then the shp file will have no info on the crowns and hence will create an empty gpd Dataframe.
- # this causes an error so skip creating geojson. The training code will also ignore png so no problem.
- try:
- filename = out_path_root.with_suffix(out_path_root.suffix + ".geojson")
- moved_scaled = overlapping_crowns.set_geometry(moved_scaled)
- moved_scaled.to_file(
- driver="GeoJSON",
- filename=filename,
- )
- with open(filename, "r") as f:
- shp = json.load(f)
- shp.update(impath)
- with open(filename, "w") as f:
- json.dump(shp, f)
- except ValueError:
- print("Cannot write empty DataFrame to file.")
- continue
- # Repeat and want to save crowns before being moved as overlap with lidar data to get the heights
- # can try clean up the code here as lots of reprojecting and resaving but just going to get to
- # work for now
- out_geo_file = out_path_root.parts[-1] + "_geo"
- out_path_geo = out_path / Path(out_geo_file)
- try:
- filename_unmoved = out_path_geo.with_suffix(out_path_geo.suffix + ".geojson")
- overlapping_crowns.to_file(
- driver="GeoJSON",
- filename=filename_unmoved,
- )
- with open(filename_unmoved, "r") as f:
- shp = json.load(f)
- shp.update(impath)
- with open(filename_unmoved, "w") as f:
- json.dump(shp, f)
- except ValueError:
- print("Cannot write empty DataFrame to file.")
- continue
-
- print("Tiling complete")
+ Tiles up large rasters into managable tiles for training and prediction. If crowns are not supplied the function
+ will tile up the entire landscape for prediction. If crowns are supplied the function will tile these with the image
+ and skip tiles without a minimum coverage of crowns. The 'threshold' can be varied to ensure a good coverage of
+ crowns across a traing tile. Tiles that do not have sufficient coverage are skipped.
+
+ Args:
+ img_path: Path to the orthomosaic
+ out_dir: Output directory
+ buffer: Overlapping buffer of tiles in meters (UTM)
+ tile_width: Tile width in meters
+ tile_height: Tile height in meters
+ crowns: Crown polygons as a geopandas dataframe
+ threshold: Min proportion of the tile covered by crowns to be accepted {0,1}
+ nan_theshold: Max proportion of tile covered by nans
+ dtype_bool: Flag to edit dtype to prevent black tiles
+
+ Returns:
+ None
+ """
+ out_path = Path(out_dir)
+ os.makedirs(out_path, exist_ok=True)
+ tilename = Path(img_path).stem
+ with rasterio.open(img_path) as data:
+ crs = data.crs.to_epsg() # Update CRS handling to avoid deprecated syntax
+
+ tile_args = [
+ (img_path, out_dir, buffer, tile_width, tile_height, dtype_bool, minx, miny, crs, tilename, crowns,
+ threshold, nan_threshold, mode, class_column)
+ for minx in np.arange(ceil(data.bounds[0]) + buffer, data.bounds[2] - tile_width - buffer, tile_width, int)
+ for miny in np.arange(ceil(data.bounds[1]) + buffer, data.bounds[3] - tile_height - buffer, tile_height,
+ int)
+ ]
+
+ with concurrent.futures.ProcessPoolExecutor() as executor: # Use ProcessPoolExecutor here
+ list(executor.map(process_tile_train_helper, tile_args))
+
+ logger.info("Tiling complete")
def image_details(fileroot):
@@ -429,32 +469,41 @@ def is_overlapping_box(test_boxes_array, train_box):
return False
-def record_data(crowns,
- out_dir,
- column='status'):
- """Function that will record a list of classes into a file that can be readed during training.
+def record_classes(crowns: gpd.GeoDataFrame, out_dir: str, column: str = 'status', save_format: str = 'json'):
+ """Function that records a list of classes into a file that can be read during training.
Args:
crowns: gpd dataframe with the crowns
out_dir: directory to save the file
column: column name to get the classes from
+ save_format: format to save the file ('json' or 'pickle')
Returns:
None
"""
-
+ # Extract unique class names from the specified column
list_of_classes = crowns[column].unique().tolist()
- print("**The list of classes are:**")
- print(list_of_classes)
- print("**The list has been saved to the out_dir**")
+ # Sort the list of classes in alphabetical order
+ list_of_classes.sort()
- # Write it into file "classes.txt"
- out_tif = out_dir + 'classes.txt'
- f = open(out_tif, "w")
- for i in list_of_classes:
- f.write("%s\n" % i)
- f.close()
+ # Create a dictionary for class-to-index mapping
+ class_to_idx = {class_name: idx for idx, class_name in enumerate(list_of_classes)}
+
+ # Save the class-to-index mapping to disk
+ out_path = Path(out_dir)
+ os.makedirs(out_path, exist_ok=True)
+
+ if save_format == 'json':
+ with open(out_path / 'class_to_idx.json', 'w') as f:
+ json.dump(class_to_idx, f)
+ elif save_format == 'pickle':
+ with open(out_path / 'class_to_idx.pkl', 'wb') as f:
+ pickle.dump(class_to_idx, f)
+ else:
+ raise ValueError("Unsupported save format. Use 'json' or 'pickle'.")
+
+ print(f"Classes saved as {save_format} file: {class_to_idx}")
def to_traintest_folders( # noqa: C901
@@ -491,7 +540,8 @@ def to_traintest_folders( # noqa: C901
Path(out_dir / "train").mkdir(parents=True, exist_ok=True)
Path(out_dir / "test").mkdir(parents=True, exist_ok=True)
- file_names = tiles_dir.glob("*.png")
+ # file_names = tiles_dir.glob("*.png")
+ file_names = tiles_dir.glob("*.geojson")
file_roots = [item.stem for item in file_names]
num = list(range(0, len(file_roots)))
@@ -540,29 +590,56 @@ def to_traintest_folders( # noqa: C901
if __name__ == "__main__":
- # Right let"s test this first with Sepilok 10cm resolution, then I need to try it with 50cm resolution.
- img_path = "/content/drive/Shareddrives/detectreeRGB/benchmark/Ortho2015_benchmark/P4_Ortho_2015.tif"
- crown_path = "gdrive/MyDrive/JamesHirst/NY/Buffalo/Buffalo_raw_data/all_crowns.shp"
- out_dir = "./"
- # Read in the tiff file
- # data = img_data.open(img_path)
- # Read in crowns
- data = rasterio.open(img_path)
+ # Define paths to the input data
+ img_path = "/path/to/your/orthomosaic.tif" # Path to your input orthomosaic file
+ crown_path = "/path/to/your/crown_shapefile.shp" # Path to the shapefile containing crowns
+ out_dir = "/path/to/output/directory" # Directory where you want to save the tiled output
+
+ # Optional parameters for tiling and processing
+ buffer = 30 # Overlap between tiles (in meters)
+ tile_width = 200 # Tile width (in meters)
+ tile_height = 200 # Tile height (in meters)
+ nan_threshold = 0.1 # Max proportion of tile that can be NaN before it's discarded
+ threshold = 0.5 # Minimum crown coverage per tile for it to be kept (0-1)
+ dtype_bool = False # Change dtype to uint8 to avoid black tiles
+ mode = "rgb" # Use 'rgb' for regular 3-channel imagery, 'ms' for multispectral
+ class_column = "species" # Column in the crowns file to use as the class label
+
+ # Read in the crowns
crowns = gpd.read_file(crown_path)
- print(
- "shape =",
- data.shape,
- ",",
- data.bounds,
- "and number of bands =",
- data.count,
- ", crs =",
- data.crs,
+
+ # Record the classes and save the class mapping
+ record_classes(
+ crowns=crowns, # Geopandas dataframe with crowns
+ out_dir=out_dir, # Output directory to save class mapping
+ column=class_column, # Column used for classes
+ save_format='json' # Choose between 'json' or 'pickle'
)
- buffer = 20
- tile_width = 200
- tile_height = 200
+ # Perform the tiling, ensuring the selected class column is used
+ tile_data(
+ img_path=img_path,
+ out_dir=out_dir,
+ buffer=buffer,
+ tile_width=tile_width,
+ tile_height=tile_height,
+ crowns=crowns,
+ threshold=threshold,
+ nan_threshold=nan_threshold,
+ dtype_bool=dtype_bool,
+ mode=mode,
+ class_column=class_column # Use the selected class column (e.g., 'species', 'status')
+ )
+
+ # Split the data into training and validation sets (optional)
+ # This can be used for train/test folder creation based on the generated tiles
+ to_traintest_folders(
+ tiles_folder=out_dir, # Directory where tiles are saved
+ out_folder="/path/to/final/data/output", # Final directory for train/test data
+ test_frac=0.2, # Fraction of data to be used for testing
+ folds=5, # Number of folds (optional, can be set to 1 for no fold splitting)
+ strict=True, # Ensure no overlap between train/test tiles
+ seed=42 # Set seed for reproducibility
+ )
- tile_data_train(data, out_dir, buffer, tile_width, tile_height, crowns)
- to_traintest_folders(folds=5)
+ logger.info("Tiling process completed successfully!")
diff --git a/detectree2/tests/test_preprocessing.py b/detectree2/tests/test_preprocessing.py
index 6f2eb9cf..49a7d383 100644
--- a/detectree2/tests/test_preprocessing.py
+++ b/detectree2/tests/test_preprocessing.py
@@ -62,11 +62,11 @@ def test_tiling(self):
tile_height = 40
threshold = 0.2
- from detectree2.preprocessing.tiling import tile_data_train
+ from detectree2.preprocessing.tiling import tile_data
out_dir = os.path.join(out_dir, "tiles")
- tile_data_train(data, out_dir, buffer, tile_width, tile_height, crowns, threshold)
+ tile_data(img_path, out_dir, buffer, tile_width, tile_height, crowns, threshold)
# TODO: install pytest-depends to automatically order
diff --git a/docs/source/index.rst b/docs/source/index.rst
index 255eefc3..545790ab 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -42,6 +42,7 @@ Accurate delineation of individual tree crowns in tropical forests from aerial R
installation
tutorial
+ tutorial_multi
contributing
using-git
.. _notebooks/contributing_guide
diff --git a/docs/source/tutorial.rst b/docs/source/tutorial.rst
index 1af9c89a..369d71ec 100644
--- a/docs/source/tutorial.rst
+++ b/docs/source/tutorial.rst
@@ -1,10 +1,12 @@
Tutorial
========
+
This tutorial goes through the steps of single class (tree) detection and
-delineation. A guide to multiclass prediction (e.g. species mapping,
-disease mapping) is coming soon. Example data that can be used in
-this tutorial is available `here `_.
+delineation from RGB and multispectral data. A guide to multiclass prediction
+(e.g. species mapping, disease mapping) is coming soon. Example data that can
+be used in this tutorial is available
+`here `_.
The key steps are:
@@ -40,9 +42,15 @@ If you would just like to make predictions on an orthomosaic with a pre-trained
model from the ``model_garden``, skip to part 4 (Generating landscape
predictions).
+The data preparation and training process for both RGB and multispectral data
+is presented here. The process is similar for both data types but there are
+some key differences that are highlighted. Training a single model on both RGB
+and multispectral data at the same time is not currently supported. Stick to
+one data type per model (or stack the RGB bands with the multispectral bands
+and treat as in the case of multispectral data).
-Preparing data
---------------
+Preparing data (RGB and multispectral)
+--------------------------------------
An example of the recommended file structure when training a new model is as follows:
@@ -58,6 +66,8 @@ An example of the recommended file structure when training a new model is as fol
├── rgb
│ ├── Paracou_RGB_2016_10cm.tif (RGB orthomosaic in local UTM CRS)
│ └── Paracou_RGB_2019.tif (RGB orthomosaic in local UTM CRS)
+ ├── ms
+ │ └── Paracou_MS_2016.tif (Multispectral orthomosaic in local UTM CRS)
└── crowns
└── UpdatedCrowns8.gpkg (Crown polygons readable by geopandas e.g. Geopackage, shapefile)
@@ -65,11 +75,16 @@ Here we have two sites available to train on (Danum and Paracou). Several site d
included in the training and testing phase (but only a single site directory is required).
If available, several RGB orthomosaics can be included in a single site directory (see e.g ``Paracou -> RGB``).
+For Paracou, we also have a multispectral scan available (5-bands). For this data, the ``mode`` parameter in the
+``tile_data`` function should be set to ``"ms"``. This calls a different routine for tiling the data that retains the
+``.tif`` format instead of converting to ``.png`` as in the case of ``rgb``. This comes at a slight expense of speed
+later on but is necessary to retain all the multispectral information.
+
We call functions to from ``detectree2``'s tiling and training modules.
.. code-block:: python
- from detectree2.preprocessing.tiling import tile_data_train, to_traintest_folders
+ from detectree2.preprocessing.tiling import tile_data, to_traintest_folders
from detectree2.models.train import register_train_data, MyTrainer, setup_cfg
import rasterio
import geopandas as gpd
@@ -83,9 +98,6 @@ Set up the paths to the orthomosaic and corresponding manual crown data.
img_path = site_path + "/rgb/2016/Paracou_RGB_2016_10cm.tif"
crown_path = site_path + "/crowns/220619_AllSpLabelled.gpkg"
- # Read in the tiff file
- data = rasterio.open(img_path)
-
# Read in crowns (then filter by an attribute if required)
crowns = gpd.read_file(crown_path)
crowns = crowns.to_crs(data.crs.data) # making sure CRS match
@@ -98,6 +110,8 @@ The tile size will depend on:
* Available computational resources.
* The detail required on the crown outline.
* If using a pre-trained model, the tile size used in training should roughly match the tile size of predictions.
+* The ``mode`` depends on whether you are tiling 3-band RGB (``mode="rgb"``) data of multispectral data of 4 or more
+bands (``mode="ms"``).
.. code-block:: python
@@ -112,24 +126,28 @@ The tile size will depend on:
The total tile size here is 100 m x 100 m (a 40 m x 40 m core area with a surrounding 30 m buffer that overlaps with
surrounding tiles). Including a buffer is recommended as it allows for tiles that include more training crowns.
-Next we tile the data. The ``tile_data_train`` function will only retain tiles that contain more than the given
-``threshold`` coverage of training data (here 60%). This helps to reduce the chance that the network is trained with
-tiles that contain a large number of unlabelled crowns (which would reduce its sensitivity).
+Next we tile the data. The ``tile_data`` function, when ``crowns`` is supplied, will only retain tiles that contain more
+than the given ``threshold`` coverage of training data (here 60%). This helps to reduce the chance that the network is
+trained with tiles that contain a large number of unlabelled crowns (which would reduce its sensitivity). This value
+should be adjusted depending on the density of crowns in the landscape (e.g. 10% may be more appropriate for savannah
+type systems or urban environments).
.. code-block:: python
- tile_data_train(data, out_dir, buffer, tile_width, tile_height, crowns, threshold)
+ tile_data(img_path, out_dir, buffer, tile_width, tile_height, crowns, threshold, mode="rgb")
.. warning::
- If tiles are outputing as blank images set ``dtype_bool = True`` in the ``tile_data_train`` function. This is a bug
- and we are working on fixing it.
+ If tiles are outputing as blank images set ``dtype_bool = True`` in the ``tile_data`` function. This is a bug
+ and we are working on fixing it. Supplying crown polygons will cause the function to tile for
+ training (as opposed to landscape prediction which is described below).
.. note::
- You will want to relax the ``threshold`` value if your trees are sparsely distributed across your landscape.
- Remember, ``detectree2`` was initially designed for dense, closed canopy forests so some of the default assumptions
- will reflect that.
+ You will want to relax the ``threshold`` value if your trees are sparsely distributed across your landscape or if you
+ want to include non-forest areas (e.g. river, roads). Remember, ``detectree2`` was initially designed for dense,
+ closed canopy forests so some of the default assumptions will reflect that and parameters will need to be adjusted
+ for different systems.
-Send geojsons to train folder (with sub-folders for k-fold cross validation) and test folder.
+Send geojsons to train folder (with sub-folders for k-fold cross validation) and a test folder.
.. code-block:: python
@@ -141,7 +159,7 @@ Send geojsons to train folder (with sub-folders for k-fold cross validation) and
that have any overlap with test tiles (including the buffers), ensuring strict spatial separation of the test data.
However, this can remove a significant proportion of the data available to train on so if validation accuracy is a
sufficient test of model performance ``test_frac`` can be set to ``0`` or set ``strict=False`` (which allows for
- some overlap in the buffers between test and train/val tiles).
+ overlap in the buffers between test and train/val tiles).
The data has now been tiled and partitioned for model training, tuning and evaluation.
@@ -161,8 +179,9 @@ The data has now been tiled and partitioned for model training, tuning and evalu
└── test (test data folder)
-It is advisable to do a visual inspection on the tiles to ensure that the tiling has worked as expected and that crowns
-and images align. This can be done quickly with the inbuilt ``detectron2`` visualisation tools.
+It is recommended to visually inspect the tiles before training to ensure that the tiling has worked as expected and
+that crowns and images align. This can be done with the inbuilt ``detectron2`` visualisation tools. For RGB tiles
+(``.png``), the following code can be used to visualise the training data.
.. code-block:: python
@@ -199,8 +218,61 @@ and images align. This can be done quickly with the inbuilt ``detectron2`` visua
|
-Training a model
-----------------
+Alternatively, with some adaptation the ``detectron2`` visualisation tools can also be used to visualise the
+multispectral (``.tif``) tiles.
+
+.. code-block:: python
+
+ import rasterio
+ from detectron2.utils.visualizer import Visualizer
+ from detectree2.models.train import combine_dicts
+ from detectron2.data import DatasetCatalog, MetadataCatalog
+ from PIL import Image
+ import numpy as np
+ import cv2
+ import matplotlib.pyplot as plt
+ from IPython.display import display
+
+ val_fold = 1
+ name = "Paracou"
+ tiles = "/tilesMS_" + appends + "/train"
+ train_location = "/content/drive/MyDrive/WORK/detectree2/data/" + name + tiles
+ dataset_dicts = combine_dicts(train_location, val_fold)
+ trees_metadata = MetadataCatalog.get(name + "_train")
+
+ # Function to normalize and convert multi-band image to RGB if needed
+ def prepare_image_for_visualization(image):
+ if image.shape[2] == 3:
+ # If the image has 3 bands, assume it's RGB
+ image = np.stack([
+ cv2.normalize(image[:, :, i], None, 0, 255, cv2.NORM_MINMAX)
+ for i in range(3)
+ ], axis=-1).astype(np.uint8)
+ else:
+ # If the image has more than 3 bands, choose the first 3 for visualization
+ image = image[:, :, :3] # Or select specific bands
+ image = np.stack([
+ cv2.normalize(image[:, :, i], None, 0, 255, cv2.NORM_MINMAX)
+ for i in range(3)
+ ], axis=-1).astype(np.uint8)
+
+ return image
+
+ # Visualize each image in the dataset
+ for d in dataset_dicts:
+ with rasterio.open(d["file_name"]) as src:
+ img = src.read() # Read all bands
+ img = np.transpose(img, (1, 2, 0)) # Convert to HWC format
+ img = prepare_image_for_visualization(img) # Normalize and prepare for visualization
+
+ visualizer = Visualizer(img[:, :, ::-1]*10, metadata=trees_metadata, scale=0.5)
+ out = visualizer.draw_dataset_dict(d)
+ image = out.get_image()[:, :, ::-1]
+ display(Image.fromarray(image))
+
+
+Training a model (RGB)
+----------------------
Before training can commence, it is necessary to register the training data. It is possible to set a validation fold for
model evaluation (which can be helpful for tuning models). The validation fold can be changed over different training
@@ -231,7 +303,7 @@ datasets should be tuples containing strings. If just a single site is being use
trains = ("Paracou_train", "Danum_train", "SepilokEast_train", "SepilokWest_train") # Registered train data
tests = ("Paracou_val", "Danum_val", "SepilokEast_val", "SepilokWest_val") # Registered validation data
- out_dir = "/content/drive/Shareddrives/detectree2/220809_train_outputs"
+ out_dir = "/content/drive/Shareddrives/detectree2/240809_train_outputs"
cfg = setup_cfg(base_model, trains, tests, workers = 4, eval_period=100, max_iter=3000, out_dir=out_dir) # update_model arg can be used to load in trained model
@@ -257,7 +329,7 @@ Then set up the configurations as before but with the trained model also supplie
trains = ("Paracou_train", "Danum_train", "SepilokEast_train", "SepilokWest_train") # Registered train data
tests = ("Paracou_val", "Danum_val", "SepilokEast_val", "SepilokWest_val") # Registered validation data
- out_dir = "/content/drive/Shareddrives/detectree2/220809_train_outputs"
+ out_dir = "/content/drive/Shareddrives/detectree2/240809_train_outputs"
cfg = setup_cfg(base_model, trains, tests, trained_model, workers = 4, eval_period=100, max_iter=3000, out_dir=out_dir) # update_model arg used to load in trained model
@@ -267,7 +339,11 @@ Then set up the configurations as before but with the trained model also supplie
model training will converge given the particularities of the data supplied and computational resources available.
Once we are all set up, we can get commence model training. Training will continue until a specified number of
-iterations (``max_iter``) or until model performance is no longer improving ("early stopping" via ``patience``).
+iterations (``max_iter``) or until model performance is no longer improving ("early stopping" via ``patience``). The
+``patience`` parameter sets the number of training epochs to wait for an improvement in validation accuracy before
+stopping training. This is useful for preventing overfitting and saving time. Each time an improved model is found it is
+saved to the output directory.
+
Training outputs, including model weights and training metrics, will be stored in ``out_dir``.
.. code-block::
@@ -281,12 +357,306 @@ Training outputs, including model weights and training metrics, will be stored i
Early stopping is implemented and will be triggered by a sustained failure to improve on the performance of
predictions on the validation fold. This is measured as the AP50 score of the validation predictions.
+Training a model (multispectral)
+--------------------------------
+
+The process for training a multispectral model is similar to that for RGB data but there are some key steps that are
+different. Data will be read from ``.tif`` files of 4 or more bands instead of the 3-band ``.png`` files.
+
+Data should be registered as before:
+
+.. code-block:: python
+
+ from detectree2.models.train import register_train_data, remove_registered_data
+ val_fold = 5
+ appends = "40_30_0.6"
+ site_path = "/content/drive/SharedDrive/detectree2/data/Paracou"
+ train_location = site_path + "/tilesMS_" + appends + "/train/"
+ register_train_data(train_location, "ParacouMS", val_fold)
+
+The number of bands can be checked with rasterio:
+
+.. code-block:: python
+
+ import rasterio
+ import os
+ import glob
+
+ # Read in geotif and assess mean and sd for each band
+ #site_path = "/content/drive/MyDrive/WORK/detectree2/data/Paracou"
+ folder_path = site_path + "/tilesMS_" + appends + "/"
+
+ # Select path of first .tif file
+ img_paths = glob.glob(folder_path + "*.tif")
+ img_path = img_paths[0]
+
+ # Open the raster file
+ with rasterio.open(img_path) as dataset:
+ # Get the number of bands
+ num_bands = dataset.count
+
+ # Print the number of bands
+ print(f'The raster has {num_bands} bands.')
+
+
+Due to the additional bands, we must modify the weights of the first convolutional layer (conv1) to accommodate a
+different number of input channels. This is done with the ``modify_conv1_weights`` function. The extension of the
+``cfg.MODEL.PIXEL_MEAN`` and ``cfg.MODEL.PIXEL_STD`` lists to include the additional bands happens within the
+``setup_cfg`` function when ``num_bands`` is set to a value greater than 3. ``imgmode`` should be set to ``"ms"`` to
+ensure the correct training routines are called.
+
+.. code-block:: python
+
+ from datetime import date
+ from detectron2.modeling import build_model
+ import torch.nn as nn
+ import torch.nn.init as init
+ from detectron2.modeling.roi_heads.fast_rcnn import FastRCNNOutputLayers
+ import numpy as np
+ from detectree2.models.train import modify_conv1_weights, MyTrainer, setup_cfg
+
+ # Good idea to keep track of the date if producing multiple models
+ today = date.today()
+ today = today.strftime("%y%m%d")
+
+ names = ["ParacouMS",]
+
+ trains = (names[0] + "_train",)
+ tests = (names[0] + "_val",)
+ out_dir = "/content/drive/SharedDrive/detectree2/models/" + today + "_ParacouMS"
+
+ base_model = "COCO-InstanceSegmentation/mask_rcnn_R_101_FPN_3x.yaml" # Path to the model config
+
+ # Set up the configuration
+ cfg = setup_cfg(base_model, trains, tests, workers = 2, eval_period=50,
+ base_lr = 0.0003, backbone_freeze=0, gamma = 0.9,
+ max_iter=500000, out_dir=out_dir, resize = "rand_fixed", imgmode="ms",
+ num_bands= num_bands) # update_model arg can be used to load in trained model
+
+ # Build the model
+ model = build_model(cfg)
+
+ # Adjust input layer to accept correct number of channels
+ modify_conv1_weights(model, num_input_channels=num_bands)
+
+
+With additional bands, more data is being passed through the network per image so it may be neessary to reduce the
+number of images per batch. Only do this is you a getting warnings/errors about memory usage (e.g.
+``CUDA out of memory``) as it will slow down training.
+
+.. code-block:: python
+
+ cfg.SOLVER.IMS_PER_BATCH = 1
+
+
+Training can now commence as before:
+
+.. code-block::
+
+ trainer = MyTrainer(cfg, patience = 5)
+ trainer.resume_or_load(resume=False)
+ trainer.train()
+
+
+Data augmentation
+-----------------
+
+Data augmentation is a technique used to artificially increase the size of the training dataset by applying random
+transformations to the input data. This can help improve the generalization of the model and reduce overfitting. The
+``detectron2`` library provides a range of data augmentation options that can be used during training. These include
+random flipping, scaling, rotation, and color jittering.
+
+Additionally, resizing of the input data can be applied as an augmentation technique. This can be useful when training
+a model that should be flexible with respect to tile size and resolution.
+
+By default, random rotations and flips will be performed on input images.
+
+.. code-block:: python
+
+ augmentations = [
+ T.RandomRotation(angle=[90, 90], expand=False),
+ T.RandomFlip(prob=0.4, horizontal=True, vertical=False),
+ T.RandomFlip(prob=0.4, horizontal=False, vertical=True),
+ ]
+
+If the input data is RGB, additional augmentations will be applied to adjust the brightness, contrast, saturation, and
+lighting of the images. These augmentations are only available for RGB images and will not be applied to multispectral.
+
+.. code-block:: python
+ # Additional augmentations for RGB images
+ if cfg.IMGMODE == "rgb":
+ augmentations.extend([
+ T.RandomBrightness(0.7, 1.5),
+ T.RandomLighting(0.7),
+ T.RandomContrast(0.6, 1.3),
+ T.RandomSaturation(0.8, 1.4)
+ ])
+
+There are three resizing modes for the input data (1) ``fixed``, (2) ``random``, and (3) ``rand_fixed``. This are set
+in the configuration file (``cfg``) with the `setup_cfg` function.
+
+The ``fixed`` mode will resize the input data to a images width/height of 1000 pixels. This is efficient but may not
+lead to models that transfer well across scales (e.g. if the model is to be used on a range of different resolutions).
+
+.. code-block:: python
+
+ if cfg.RESIZE == "fixed":
+ augmentations.append(T.ResizeShortestEdge([1000, 1000], 1333))
+
+The ``random`` mode will randomly resize (and resample to change the resolutions) the input data to between 0.6 and 1.4
+times the original height/width. This can help the model learn to detect objects at different scales and from images of
+different resolutions (and sensors).
+
+.. code-block:: python
+
+ elif cfg.RESIZE == "random":
+ size = None
+ for i, datas in enumerate(DatasetCatalog.get(cfg.DATASETS.TRAIN[0])):
+ location = datas['file_name']
+ try:
+ # Try to read with cv2 (for RGB images)
+ img = cv2.imread(location)
+ if img is not None:
+ size = img.shape[0]
+ else:
+ # Fall back to rasterio for multi-band images
+ with rasterio.open(location) as src:
+ size = src.height # Assuming square images
+ except Exception as e:
+ # Handle any errors that occur during loading
+ print(f"Error loading image {location}: {e}")
+ continue
+ break
+
+ if size:
+ print("ADD RANDOM RESIZE WITH SIZE = ", size)
+ augmentations.append(T.ResizeScale(0.6, 1.4, size, size))
+
+The ``rand_fixed`` mode constrains the random resizing to a fixed pixel width/height range (regardless of the resolution
+of the input data). This can help to speed up training if the input tiles are high resolution and pushing up against
+available memory limits. It retains the benefits of random resizing but constrains the range of possible sizes.
+
+.. code-block:: python
+
+ elif cfg.RESIZE == "rand_fixed":
+ augmentations.append(T.ResizeScale(0.6, 1.4, 1000, 1000))
+
+Which resizing option is selected depends on the problem at hand. A more precise delineation can be generated if high
+resolution images are retained but this comes at the cost of increased memory usage and slower training times. If the
+model is to be used on a range of different resolutions, random resizing can help the model learn to detect objects at
+different scales.
+
+
+Post-training (check training convergence)
+------------------------------------------
+
+It is important to check that the model has converged and is not overfitting. This can be done by plotting the training
+and validation loss over time. The ``detectron2`` training routine will output a ``metrics.json`` file that can be used
+to plot the training and validation loss. The following code can be used to plot the loss:
+
+.. code-block:: python
+
+ import json
+ import matplotlib.pyplot as plt
+ from detectree2.models.train import load_json_arr
+
+ #out_dir = "/content/drive/Shareddrives/detectree2/models/230103_resize_full"
+ experiment_folder = out_dir
+
+ experiment_metrics = load_json_arr(experiment_folder + '/metrics.json')
+
+ plt.plot(
+ [x['iteration'] for x in experiment_metrics if 'validation_loss' in x],
+ [x['validation_loss'] for x in experiment_metrics if 'validation_loss' in x], label='Total Validation Loss', color='red')
+ plt.plot(
+ [x['iteration'] for x in experiment_metrics if 'total_loss' in x],
+ [x['total_loss'] for x in experiment_metrics if 'total_loss' in x], label='Total Training Loss')
+
+ plt.legend(loc='upper right')
+ plt.title('Comparison of the training and validation loss of detectree2')
+ plt.ylabel('Total Loss')
+ plt.xlabel('Number of Iterations')
+ plt.show()
+
+.. image:: ../../report/figures/train_val_loss.png
+ :width: 400
+ :alt: Train and validation loss
+ :align: center
+
+|
+Training loss and validation loss decreased over time. As training continued, the validation loss flattened whereas the
+training loss continued to decrease. The ``patience`` mechanism prevented training from continuing after 3000 iterations
+preventing overfitting. If validation loss is substantially higher than training loss, the model may be overfitted.
+
+To understand how the segmentation performance improves through training, it is also possible to plot the AP50 score
+(see below for definition) over the iterations. This can be done with the following code:
+
+
+.. code-block:: python
+
+ plt.plot(
+ [x['iteration'] for x in experiment_metrics if 'validation_loss' in x],
+ [x['validation_loss'] for x in experiment_metrics if 'validation_loss' in x], label='Total Validation Loss', color='red')
+ plt.plot(
+ [x['iteration'] for x in experiment_metrics if 'total_loss' in x],
+ [x['total_loss'] for x in experiment_metrics if 'total_loss' in x], label='Total Training Loss')
+
+ plt.legend(loc='upper right')
+ plt.title('Comparison of the training and validation loss of detectree2')
+ plt.ylabel('Total Loss')
+ plt.xlabel('Number of Iterations')
+ plt.show()
+
+.. image:: ../../report/figures/val_AP50.png
+ :width: 400
+ :alt: AP50 score
+ :align: center
+|
+
+Performance metrics
+-------------------
+
+In instance segmentation, **AP50** refers to the **Average Precision** at an Intersection over Union (IoU) threshold of
+**50%**.
+
+- **Precision**: Precision is the ratio of correctly predicted positive objects (true positives) to all predicted
+ bjects (both true positives and false positives).
+
+ - Formula: :math:`\text{Precision} = \frac{\text{True Positives}}{\text{True Positives} + \text{False Positives}}`
+
+- **Recall**: Recall is the ratio of correctly predicted positive objects (true positives) to all actual positive
+objects in the ground truth (true positives and false negatives).
+
+ - Formula: :math:`\text{Recall} = \frac{\text{True Positives}}{\text{True Positives} + \text{False Negatives}}`
+
+- **Average Precision (AP)**: AP is a common metric used to evaluate the performance of object detection and instance
+segmentation models. It represents the precision of the model across various recall levels. In simpler terms, it is a
+combination of the model's ability to correctly detect objects and how complete those detections are.
+
+- **IoU (Intersection over Union)**: IoU measures the overlap between the predicted segmentation mask (or bounding box
+in object detection) and the ground truth mask. It is calculated as the area of overlap divided by the area of union
+between the predicted and true masks.
+
+- **AP50**: Specifically, **AP50** computes the average precision for all object classes at a threshold of **50% IoU**.
+This means that a predicted object is considered correct (a true positive) if the IoU between the predicted and ground
+truth masks is greater than or equal to 0.5 (50%). It is a relatively lenient threshold, focusing on whether the
+detected objects overlap reasonably with the ground truth, even if the boundaries aren't perfectly aligned.
+
+In summary, AP50 evaluates how well a model detects objects with a 50% overlap between the predicted and ground truth
+masks in instance segmentation tasks.
+
+.. image:: ../../report/figures/IoU_AP.png
+ :width: 400
+ :alt: IoU and AP illustration
+ :align: center
+
Evaluating model performance
----------------------------
Coming soon! See Colab notebook for example routine (``detectree2/notebooks/colab/evaluationJB.ipynb``).
+
Generating landscape predictions
--------------------------------
@@ -312,8 +682,7 @@ can discard partial the crowns predicted at the edge of tiles.
site_path = "/content/drive/Shareddrives/detectree2/data/BCI_50ha"
img_path = site_path + "/rgb/2015.06.10_07cm_ORTHO.tif"
tiles_path = site_path + "/tilespred/"
- # Read in the geotiff
- data = rasterio.open(img_path)
+
# Location of trained model
model_path = "/content/drive/Shareddrives/detectree2/models/220629_ParacouSepilokDanum_JB.pth"
@@ -321,11 +690,12 @@ can discard partial the crowns predicted at the edge of tiles.
buffer = 30
tile_width = 40
tile_height = 40
- tile_data(data, tiles_path, buffer, tile_width, tile_height, dtype_bool = True)
+ tile_data(img_path, tiles_path, buffer, tile_width, tile_height, dtype_bool = True)
.. warning::
If tiles are outputing as blank images set ``dtype_bool = True`` in the ``tile_data`` function. This is a bug
- and we are working on fixing it.
+ and we are working on fixing it. Avoid supplying crown polygons otherwise the function will run as if it is tiling
+ for training.
To download a pre-trained model from the ``model_garden`` you can run ``wget`` on the package repo
diff --git a/docs/source/tutorial_multi.rst b/docs/source/tutorial_multi.rst
index 21b1fc72..125e06ef 100644
--- a/docs/source/tutorial_multi.rst
+++ b/docs/source/tutorial_multi.rst
@@ -4,9 +4,10 @@ Tutorial (multiclass)
This tutorial goes through the steps of multiclass detection and
delineation (e.g. species mapping, disease mapping). A guide to single
class prediction is available
-`here `_. The multiclass
-process is more complicated than single class prediction as the classes need to
-be correctly encoded in the data.
+`here `_ - this covers
+more detail on the fundamentals of training and should be reviewed before this
+tutorial. The multiclassprocess is slightly more intricate than single class
+prediction as the classes need to be correctly encoded and caried throughout the pipeline.
The key steps are:
@@ -15,4 +16,161 @@ The key steps are:
3. Evaluating model performance
4. Making landscape level predictions
-THE REST OF THIS TUTORIAL IS UNDER CONSTRUCTION
+
+Preparing data (RGB and multispectral)
+--------------------------------------
+
+Data can be prepared in a similar way to the single class case but the classes
+and their order (mapping) need to be saved so that they can be accessed
+consistently across training and prediction. The classes are saved in a json
+file with the class names and their indices. The indices are used to encode
+the classes in the training.
+
+.. code-block:: python
+
+ import rasterio
+ import geopandas as gpd
+
+ # Load the data
+ base_dir = "/content/drive/MyDrive/SHARED/detectree2"
+
+ site_path = base_dir + "/data/Danum_lianas"
+
+ # Set the path to the orthomosaic and the crown shapefile
+ img_path = site_path + "/rgb/2017_50ha_Ortho_reproject.tif"
+ crown_path = site_path + "/crowns/Danum_lianas_full2017.gpkg"
+
+ # Here, we set the name of the output folder.
+ # Set tiling parameters
+ buffer = 30
+ tile_width = 40
+ tile_height = 40
+ threshold = 0.6
+ appends = str(tile_width) + "_" + str(buffer) + "_" + str(threshold)
+
+ out_dir = site_path + "/tilesClass_" + appends + "/"
+
+ # Read in the tiff file
+ data = rasterio.open(img_path)
+
+ # Read in crowns (then filter by an attribute?)
+ crowns = gpd.read_file(crown_path)
+ crowns = crowns.to_crs(data.crs.data)
+ print(crowns.head())
+
+ class_column = 'status'
+
+ # Record the classes and save the class mapping
+ record_classes(
+ crowns=crowns, # Geopandas dataframe with crowns
+ out_dir=out_dir, # Output directory to save class mapping
+ column=class_column, # Column to be used for classes
+ save_format='json' # Choose between 'json' or 'pickle'
+ )
+
+
+The class mapping has been saved in the output directory as a json file called
+``class_to_idx.json``. This file can now be accessed to encode the classes in
+training and prediction steps.
+
+To tile the data, we call the ``tile_data`` function as we did in the single
+class case except now we point to the column name of the classes.
+
+.. code-block:: python
+
+ # Tile the data
+ tile_data(
+ img_path=img_path, # Path to the orthomosaic
+ out_dir=out_dir, # Output directory to save tiles
+ buffer=buffer, # Buffer around the crowns
+ tile_width=tile_width, # Width of the tiles
+ tile_height=tile_height, # Height of the tiles
+ crowns=crowns, # Geopandas dataframe with crowns
+ threshold=threshold, # Threshold for the buffer
+ class_column=class_column, # Column to be used for classes
+ )
+
+ # Split the data into training and validation sets
+ to_traintest_folders(
+ tiles_folder=out_dir, # Directory where tiles are saved
+ out_folder=out_dir, # Final directory for train/test data
+ test_frac=0, # Fraction of data to be used for testing
+ folds=5, # Number of folds (optional, can be set to 1 for no fold splitting)
+ strict=False, # Ensure no overlap between train/test tiles
+ seed=42 # Set seed for reproducibility
+ )
+
+
+Training models
+---------------
+
+To train with multiple classes, we need to ensure that the classes are
+registered correctly in the dataset catalogue. This can be done with the class
+mapping file that was saved in the previous step. The class mapping file will
+set the classes and their indices.
+
+.. code-block:: python
+
+ from detectree2.models.train import register_train_data, remove_registered_data, setup_cfg, MyTrainer
+ from detectree2.preprocessing.tiling import load_class_mapping
+
+ # Set validation fold
+ val_fold = 5
+
+ site_path = base_dir + "/data/Danum_lianas"
+ train_dir = site_path + "/tilesClass_40_30_0.6/train"
+ class_mapping_file = site_path + "/tilesClass_40_30_0.6/" + "/class_to_idx.json"
+ data_name = "DanumLiana"
+
+ register_train_data(train_dir, data_name, val_fold=val_fold, class_mapping_file=class_mapping_file)
+
+
+Now the data is registered, should generate the configuration (`cfg`) and train
+the model. By passing the class mapping file to the configuration set up, the
+`cfg` will be register the number of classes.
+
+.. code-block:: python
+
+ from detectron2.modeling import build_model
+ from detectron2.modeling.roi_heads.fast_rcnn import FastRCNNOutputLayers
+ import numpy as np
+ from datetime import date
+
+
+ today = date.today()
+ today = today.strftime("%y%m%d")
+
+ names = [data_name,]
+
+ trains = (names[0] + "_train",)
+ tests = (names[0] + "_val",)
+ out_dir = "/content/drive/MyDrive/WORK/detectree2/models/" + today + "_Danum_lianas"
+
+ base_model = "COCO-InstanceSegmentation/mask_rcnn_R_101_FPN_3x.yaml" # Path to the model config
+
+ # When you increase the number of channels (i.e., the number of filters) in a Convolutional Neural Network (CNN), the general recommendation is to decrease the learning rate
+ lrs = [0.03, 0.003, 0.0003, 0.00003]
+
+ # Set up model configuration, using the class mapping to determine the number of classes
+ cfg = setup_cfg(
+ base_model="COCO-InstanceSegmentation/mask_rcnn_R_101_FPN_3x.yaml",
+ trains=trains,
+ tests=tests,
+ max_iter=500000,
+ eval_period=50,
+ base_lr=lrs[0],
+ out_dir=out_dir,
+ resize="rand_fixed",
+ class_mapping_file=class_mapping_file # Optional
+ )
+
+ # Train the model
+ trainer = MyTrainer(cfg, patience=5)
+ trainer.resume_or_load(resume=False)
+ trainer.train()
+
+
+Landscape predictions
+---------------------
+
+COMING SOON
\ No newline at end of file
diff --git a/report/figures/IoU_AP.png b/report/figures/IoU_AP.png
new file mode 100644
index 0000000000000000000000000000000000000000..603778b1c02958a38316db332b5748f78bc3cf8f
GIT binary patch
literal 251986
zcmeFYWn5d&_bwXTokF0v7Y*)STnZE|PH^|)?he7VEzlO%(&AbuUZl9YyIW5BwV(XY
zxi8MWFE5|4cQSi)t!JLKX3gwqHB~t*bW(Hx0Dz?+|55`0K!klo0H7knZcRmJHX3T$Xw{9q^yIoqvI=8Hc26FPF7nX
zPA*PfE*=mk7dIzAuK=$gh@X?yU5JYt#K*tGT6+#!K11vci6eGX7Qj+qZAo-*U4%xms~>3JMBxfVeog
zxY%Gh*xbAw-A%pN9NlPPHT`E*FD=}_t~So@HcpOIzpFMibMkN(Wn}y(*PmT72mjNI
zvxlqwAB~uUIV|igeoM^_HW1GL=%qPW$lb=>-r~PHE@|)nUn%xBzegfuZ|Z0z%IL*r
zZeeNaVeif;rlMitMx|ll?BwcB#YV*`!tuLv_BQ`k`j3?VtdQe>)9fGF{%w^1mjaM<
zcXzch`$q`BHTg%1$8X88|JJ9e`EN0R{zD9ZmJ*UQcl)=F{~pO}3zz>&feGUG17}lL
zHw#T~XA4nAm;{_WT)`IquKm9%lyY|;02kPasN}-pGE#hO|oXNig>wr
zx%hYmIe9qw1$hNP|CIUT;r}S2>Er>n{yi))u77g=EA3y|{|{>ZPjUS{)c+Lg|3`{;
z{oG|+M4`lwqJ^#)z{|+HoW}!duQf^W!8Pc=T|2jLx?AD;L#}25cmiDZLgvL4WE6_(X
z0sJ(6=Mqwzoz&9yNOo9U5Pv%UAdyl}I?iQo`&G}J!);4*i?y~u#`CpJiSl(>ZE38*
zs8emppt4>umzDh>8Mg6uwLe5gUu?q^?kNpcR|jK-#_W^|tTq+HVhtM-$UO)A$2V33
zi=~aV6IZnp6c#^Jdkrk?n@MUR9<75pyVC|ERj#i;%}z2)M`$+K=+Uu})D5j)S>(ry
zZ&{1<2MDVUCcTfp_JtS8)z-TFDGh5PP$#Rl-?r_?uLq&&^ZjzSEQgBPl6^Y2Ue>y4
z>^iaph^FJd$i7cO*{xRP)`YU2Rq`V1JS?4~ZaRO<#48`xZrR~0@8(A~;cl*15l-jU
zUdgRjbXF6TB(*?=%HV_gP~Q9)|NfjFNXQ1{#1`Q5irSxuErkV^pVUvvVUrK6N9Ev
z#~v(+kLX;oOA`r_2L5mq9=g4-f_3-)kV$S`A)yy?@TGf*chmh&FfBdZyi{dt;qdKi
ztHz2tSY1}@3i)x%oA3PdKdoO2k{R8p0B&}kU{ebUJ!sXdrl0o6zddZ_J|B8pE<9-v
zPC`WmpG%cJ)sfzmgDuBHUbGlv{grYt@d`pQ{~opdXCnC+$D@gyJ1!ltR^LG>*UY$j60kGv!>z(N~i6dN*!n+wJ-2K>Z9;5w`?vWPasOG7=>2@Bt
zcT6!KU>`s07e%9%xA@bsl%OOB;Qeh*SQyX*BUnP=g*5suL@7ffB8ZxsP1aSyHhXCH
zb%)t5`<;?@~L=-arRDM#GK(t69E3`0MkHZO2EvT>l6P1vr$KOSc9;tr4XWs
z|sNOWD=5`=j^w-4)oz~qD
zU)+pZb9+0)!NGy?$!vHqG4Br|l~TFP#uk+(*umA42V`u`Bb?V%4qp>gpQ{mWX2fR$
znE&hqC0Cf^%u3e|#|*P^v4xIJqm~9KDis%Ely8tWQ|^q
zE+|=nx7~6}uNEgI8#&8H=2uXU+A`pWeTSoEMx^Fz#5kA9LwJ*!bK==+C!KA(iEH7w
zF1TV$-Sz36Zv@u(66&eO45a`yLK$_~M@;jgo6jl!DI|W4o+yQzXNWTnE}7{PM?*E_
ziHbfyMIk)Em57^UTSq=|lBZvpQTvu85hXRtAHnFXpAmWqS~Q9OzQRC&f?PXjbW4<#
z1?Frd?(BP%x1RTy8aD=03&zs7zD`tBC)`B+SDQS9}BJ#X4=evgw-=2ndg73GO}T6Wn=
z&^zhN24_PoKkjL298(uA^K7dmvrN#ob(PrH)~^s1u3Kf0uGE+BO4wTPE7qR_-3Je{
zK3Nf9zQzd|atGOuUw-|Z<=u;g6!K}Czh63|vc2E+j_UlVUKJgt+Jp~ifImb_xG613
zf&DEFiYh&@dr=d?JUjvt_r=s?|EOQ0Ix{{i}vTChtkymSysuj
zT|u+<&SbZf*G;zZoM{}A3d>Hf2a+Q6iir7u5&nC83~(<21V4%g9le9|m(9flhZXNQ
z=N$cAabS2lMo|I$IlGwNiiZ5W{5msaysr7kZQd{(;2lmyAC3tSrtPDg-W*o82;7?t
zyy)5bp^euGl&)G>`MKNyR2r;InpO1(0O2EC0T%{&Zp3TZwva+PZqOwaeKv{MzJ?@+
z1vJXA<>fT6gtc&C*L$53eh-q0AREnpx1Ip*wMCLqPrZ&}0Q`xr)aHDRhy$VwCUg0T
z%Aagb8Uv_GvbgV$olN5jZ8BEn-?7TYL6_oPLOL{|M5Z08E-v`aPH9f#s@l7BjLh3m
zMRlt6qOS5TlU{-96tb{Af0&}dIW(9hIrP{A`(9{tA7x)4pp#JH697tV#$C(Qv#5~!
z!x}fdESNqs#hfkuL=5`PG)ytBI#j6M0_I=8^?(kNGlj^%$_ndQ5n|6W8+~#l0tvJb
zYF!$E%t*qeTQ7qekGGI~BzSB6H(^!4q&JDW^5Vsr!uiQoa
zlH!?+dH`wmKO~{nve7dB_{xAJAn`mNhB;A{uM4n}To=9v89I%Ecbp|9ke?d3T^G~3`Hzy+>?T!s}+uFADnzW0gB1OjWEphg;0#IB+m
zlXBGRqPTg4as|ncfeEjs*>oYC+p5$w2OGpC^GICtu<;|q=WPh&nTs-rz>pQ7`Y$75
zf*KtA=;Zbu%0^)w5FCEN>fbIy;WIfzW;0%7xqYI6+XKl&rr|0{dye@6Yij{bkpz#a
z3|&Oo_tL`bh`5BtKP&YFj$45Oh?$OAl-VM_NAG+uA;CQ>q)NEZ#@2w)vyUW#(5}R8
zIOmwxGwCR-^DvzZ6W0t2yPNL9JJSs3A_?Ve^9Gd(w#Bi6$PY~IbN^zT4&a@?kh!3W{
zR~+a7OXY`BJ^&cM9FB6pE59br#MO|Si6NbS^1>Q3uD4tl7ICp2zTiS$|fmtzPrnu+5m*7M#U{P
zA;#YWJ}=t=e|=QgVGAkQxU7sBml5HV@I+K>Ar#;$uYo*yQ%ihVuGb38X2m5zdDNmZ
z7Cg(T88wb7g-?Eq*28#utGxD(lyc{DT~a=MaRR;I8Ud)a;FN&N{9`xO!8pfP>lq~5
zhHeAb-0Z`J63Va-Oq6*4DDq((4Zj2NV5p#=c;jJXC5NnR@fC?{93haP>}iUQ;I^ex
zEv64Y3JpGi=ju5a4-NBCZ1Ie!@q*#Zi$mJ;EY^KcfAUb|&^?bzusF`T$%51oFOxN4Vk*0TRxO`x|o=Q5c)uhX2DQXG-9ex88>*9N#9PVq%zr6MiWu!cfsX_q^w!
zloM=xg0nY>VyU{Uk|Gfz+N)S}A0u#2-cy;lFmwX>vGQe@EPAeXXmSdtKe$jLqviYH
zmu(^hN_4ZY7O!lr(Btpv<~S!Un^TJpYj{$GO82ZQ5BRx8H07*Ux5#{5=vikemawM;
zVn9Ewwq|7`nt{^xI!qYnLR?ZzxJV?1?B(Fhm-y^xx!bN#!J7xGN&B2>9B7Pq`)#}H$6Gqz+h<%8
zy(zMVZ8iLqL186YVjV3~*$X7#K+*)uSLdSo6tLVR5}HHHarRHSSAXHlWxB~lc(++0Q6T78)ZOj8Qx8szj(wR
zHb;y7?Ht@z75+<6wcJfTM?0hihPPeHCzORZTyZlbt-}UBA&nb!9?d}M^_p%mEo%`+
zl+Om&$u;;v0!Tk~BQ~iET|4|9H4mC`;JCTBMw7=?r<^An=K>R2$2|5tgD3fL(c{o`
z!JGn^j^i1W7Pp>nHmiEiNTTmOQ`X&6JbFvRY?`{x_l^qO>~4g_2zB?Z8!gi9w&m^P
zY1c8R6DZ4J-xui|B1mFqxTQB
zHWnTvuybaVk0~a8uk=R?GKr_hL4YCZhcYB`8w%&&9vTLm(dVv=PRzS1Bovn&gu^oh
zUd;|ysohZ*$H}YYc{XYH-xWTH4;(pc6051bbx`Fkkig9?(Ms{{G+S%$mE{lQPWXWE
zp6l_s@EwkHM(Y?asbj+J0V14n1usQGq?q{9Q{u&=*tNTuT=y_H({&@~4X1v0?1DH{
zpU|}9CnbfyE>{~kmOEYBhm`?GK8sf->e*A)HjRnH0Mo_Jz33ZD603_>zAz)pcNd8E
zmm`~yGE0!ZOxhBWuzLoatPtsit9wA0PL*`0ip|YIrTjDRE%vHP{tP
zjAHzPCEM0lb0_2>EcfEJ$8_1Gj6K9f##mmZjOz#^3g%}3*PR()-cX_H9oOkrA5iJ*
z!kXMXG}+WUifVx$`aJj0FB?|`F0uP;V5ESXsoBqGvtV4Fh|xBDbI#jn>Cp`)SHZ=A
zkeo*yOF{&#nD;5Sv8cQga8@LYBE8qRFOxudL?$XyjCzHR#K81>nExg$D7!~0w}20G
zAuPwTsDC?#AD9Nn>A$G3?iN8GZLg8qJPrfbr&c^*@DlDGj?i
zBsnLo{wHtcgI-lWP}s#se)Bt#mbFL>^$Rui?)Y_det=zdcG49nw>Z2m9gQ=_RGrLP
zV<0@bbiz!XYTjU2@xli#3Q;PlpL8wT|FU&24p{d5rCClWkZZnZYm50;h+3FfSW0-&A$h9WP#Vf}q%scbaP-nER_2LH(D+1raJft_UwIY(NX2Upz+w
z+%i+asx{{=v4I+&oFCR*OF9(dQ)aUGQ$zd
zc>l{z9Wlw0g|30LJIU)I&om~?z|O(pa|#Lk6x~^blztZ9bQ-T;g)9S;5(f|}J&nzG
zngw#G*U;W?SRQ+*Hc--X>iyb(rKJM3!HBVn^MBHOyyCTOP%*2^LVT%=+d~ur3_=O>
z-F^}qq4PSu81%XLjqZUPHhvDDb88(k&*+nrt_?&U}&~n$z{1Su5d3)}g8nX_^pN
z@L?h1pB0j@qM04$16=o6l4Y#~i=kO#K*I2oK2OT`XyOiF3C;Pe#FEm;m!m{HwA>uu
zl%{=Lsm`k%eBo}TqL01NS73tkwi0)7?V(SO@I9U(d*RgYx8^1dd^H}$Yn_yzk*Onv
zsjc<9@PWmI`w=h-w%za7gI2xn)rs~$*))AuIB_abF0?3M!+<~2#Z95io`M11iYf2
z*;McsU+0wX+5k>OKVD|%19}kcWkIe=DQI~`S<%*3YqVWi7FjjqDrxZfq7d|wGB5
z{pbT~xFIGgq+bUoah+SlIOscjc9Fs@ly#
zw3()7zd`KK7u*l#i%Oq|ZcCruexADqFZ+_9(*bk)8E9`iR_x=`w*1p@0KgR6$rXW>
z0u|=C(7=w&{`z3k+3(q5tTx|E8%vl^xE$m@i@-6)wW5EJH2$1>LoNwVB%P$#;6$Z)
zjf~6@43j?p^6KSc`*#u)GhN5r;-fjt0slP(6g!Q$;HEt3@}Y7{_mvNlL%ZVvq&6Z*eAqKu9=`_x8DcPx9J`95=5=;v6R
zGg$d9QEXgxKZ&J(^-M})KCT6MHo;4n*q>QhnL|y-kvob~0(j6Q=o7!vL5d_^AdPJ6
zZiCDcknel7*HJ87D?bB2=fahlOQTXNp3rM*-Op@`i)C~ttAVM7k(0nF@vWvasJh&Om|zi<;Q8t$+rg~;qKmSlq_1zJ~Gw{;!ycv1qP>xApD$!Btg__zYE;#A;;R
zfS4E|kdpw?a{4yeb|Ei_6Rm9co~k45k<#UHg>ga1FZRfq5qX#Up5hk0hxy@~7e~Zt
zK}%QtAZ3j_>MvpL8YsEXX%zUUAV1o14ar}b@4fkehjjTPx&
zd6OMq`^(CeRz^;kUlIIeD-<@3@L$JLz~fheopslL*h64xto$R~OxQb#SE!wkd%R-G
zkW(NhW#-Lh
zzHx8q!wlVM(DxedmJ$8tWyWiNv4*qdlo?X**-u$)n&`oL5b;uqhn*=W7eqfTm#xU0_f{qfxeRAvgwd608T@cV#Vs(jfd)b9bsBGtP4&J5IX`Us!{94BYc+)Fy=
z^%X|UJrwiG;4f$mXEI59w@9ai0rKrLh3WObFxAuv7bqV(VIP!M)V~NrC$YzJXnOiv
z?KO>Xs2yw^(G@DN?mw4FY^+er%OpL&T!+b^db}m0c*_P~zdIX9#@L{0#+~aSC%d;8
zAQnnpLo)0hEXyTlYX|ee4TRci`Bvwi22(Yh16C)>zZ`=Dy+X@Ss+?QSP0|E_zIe)p
zcUYW_CMsMu2*5;q%=fvCF~#Zb)xEWg?NYRP;1ZWoMf{hcou}Tt6L6SBV6qW&x&~QW
zJHP=B|ICr`2UOa5&l6?j{(w2DKwa9_^53~8k-2pq0VTq&X{rc4g*Gt1R(fDbW~~0;r+^mNEJSTAL<cI=L9516R^$t`o4yw>9gqm%=GdhmrNWjU@wi}b*W|EdPW9=n*316H#
z*1!CX`Dex+p^Jl&Vn{0F9$1DLQ|biOJ!iG6@I8GWw-wcvONy>K_q^0S95IB%7Xeig
zKFh(?>%Y~7Yek^;Ea-!uc_SDt2!ojYq{7L6
z-%#9lujJF;+8(=^-h5zh^Aa<-;V+;BA?oECW0O-gqG_#dJv
z!J!>NIP^qHxFsnVL~xJ(jupH!m&xOn&(D$fjnc!-IZV
zw+pb3#a7J$&5)L3`$*)>Y*{B?)&6Rg;#mIldOw+-&EQu|5MoQJfC3DQlLlqk8t<3}
zZ``4_`>(MQm&aQMGZchN^R#P#Fq+
zv6!2MZ$S5ncY!nj>{9xhED@Zz*H^^R4z1#lXo>h0@HP85UpH1UzyZO1U-@=DqG*x9
z)f{W5ZSOEGa_vF>z%X3kGf~;LM9%afUewN4gm6oAy06V#bj>=9Z$uc#4G`=MHtBJv
zK}^_s)?xjy1)7-8LNfGDCalW;&`$jFJcr+@6z5Ohe?TrQw3DdSuYfz2GVxqNZoyG&r_f{jXVHzq=;}cN*~3^#pPa*
z51T{dM)KJWW8D3^%Cnt5qet6_jst6PFvA48WqW3M4Jot0IzNblQ()1lW#pzrv-L@g
z5b-puR7iYREHlfBNsqbS94RS#G0;s|8-JIDZ>f@=z%?6Q(pC&A+)cDuNt3GJN!x&rh{lPx
z5(^wkgkb-=RwD{Kg0pvkI8bo4lUqAr%R3#WI8r1XTvdm{`vr|dkmYE1tEb!hE4pJ3
zu=rr=Ov06eYCR#i@Vp;@KHDrDBDaC_`nu9KvFq&RCEyKZ3_N4Bx3{luYVu)ZFZmmg
zlM{T4W9&L_jZM|a#)r1_?91hDTY!o$N=v_GT|6eToJ=4Fe+P7b7_2@}S9(*!Vf?Yy
z56n7ZIjN$UB+(hPfK&EbXlvh&m1nA~ga3%q@E$TS){Y**7=HR0Icn1@G*jU=|8W&L
z^n;Sw>-0B_g*gExW_&x+(RGzvNIYH-biq_{+t|vRP;i(8NsKM#JL(%U&$l)D&EGnc
z>X}Mj%i#`m8|I2qmHcAb)(?9V^8%q)U}_OP=Y7#vUr1Z~h~+7hQU&&)UO{oFN;e|(8cmEt96Hba
zJy8R1IiQ>Otwn$B-w;m>5P@`!nlBCUGh7^Z*8w+bH2A>u$AJZQd)8j*9TiQ-QC+CU7rTvjU&a_k+Ac++tTOF
z!JQ9>SvsWU7pz~#3PO+F&gm+t5(N6nx4W<~BB_<%mpESd7w$QTVl(xpPK2e8$}OsK|93+J`%Q5n)lll?tSS8LQm*
zk=v{aQ{@BbE5mnrZ0h>f!NES`G#*&FDrP6*>bSs!?A5{M7M|3$mCgM?=@&?h+mtJJ
zPx~$Im-?_p*O4JO&>4dgee1bhu908ZuMt5-W3rOdzJN!ycs5Yj_)lR>vIDgF`G_2S
zlE}z(E}Yaz6S)N0+95N$AG_h?s$n@{@$6O02CKu}z5*A6m~^yCSS;tX6oQ_ykz18(
z>3c1$Z_|(4z76aX3vYC;##KMe%E7j`LciExO!BO{b=k+_phDB7nZc
z?m5!(N;B1h2+Q+~D_>Zs2c`R770{CN039@b;-b8t>GOIx?eks@CuW#0@W0Z_4}9{@
z_gag-_SO5#qpOtuRs&@iHGT@To=$JPo?!%TpkT(Z8!7rRNp`unP*69mphPivm}*d7
zYrva0Xx#4QMBAE_FI>NegiePXm1z
zX-rr{ACmHxJOU0UkIy>hAWlccHRL5#J*wLY$-RB(|5DuP7>OySsxqy(2
zt7&&W;X^Gp4A(6KC9OXArR
z=a{6p>=Fe|oIZT-$bAS_5GQ2dQ6P$+8`-zh7d>L8$1J|vmr`yMsniI{({!uz>l*mr
zR)1{#(eQR%B$x__ho9pzTJx%w?=22Zm+CPx@ASsSS|JRFncDu*lOGM1yu*(TZxZdv
zjV%11@)ePAPw<)Cx;ra3H;LVh?$k}LEmidprGk8D$efU@(Tn*nFuI>Q#s=9IgEvt!
zbaojM)s^f!#rv}-JasBb4H;~OI}fH==FBN6BvTSlaPm$$W!8&CtAy0m$5@t@zUe@_qJR
z4AhF$9#TFQ2EI(`nRr?cd`P_2f$kZT-#$EFzCng=Ek8UT=kXWM(EigP9?_r|oS%JZ
zH`T~&g
z63=ErVnypEBM)SZBS7U6j1+WL;npkQOrlp2^IYRsFpQ_lzrABF%R67FJ44JW!=w8-
zdJWt#q(hw&2vKfDjKN!IVfDMO!dSXfQ0*qnHNN7}7zLpmgW+j7u*4*qEq5m2!y^}-
zW6LEBE^^zG)5V3XKV3y>6ICPUK1&$g^up1(;!e-Whnx#WKYRJ%YlmXl3S$mdz`~mi
z=Sr1W1We8`DfJ+?OwZ@%Ghq2o2pnxuBG+VF{fx;Ug35GOzf7+Lqksx-G0V#rqK%1R
zz7W+?Mt9$I`_tuH`+FGnc)TA^x^dbI=K|MzvGTpQvaa3!p3>~=TTl8hAZ)(B9IgFn
z0Iqr)o7cUM8ABnHMIe<)a?URd@B9$Zn;f`Sc?=RjzNLPCpF#306j-j)7e7iIqLxZ~
zW)_6gM%t(2jbXCB5tL;&u#ufM4;7xWU?>D&MNIx;sG>9oO!3tb3;y1ljvR*Zoa1aR
zz8&mvL?Ltuh9HJL#`YnOs;lR-OtRRMn2&C8p`8Ad|2
z4%$!RTQEvI6Hf4ltpS$yEnxfqFp#A-$F-21?0tcDnYrabFs|$7+rJwHLevE_39n_o
zIc$j@-=t(skGz~{`#|R6anHNdbV=sq)H>8@LR`Ay=ZaZlB_0+OrrY|KFTU`Oq3|ie
z79-(;aN#4|w?$kDhG+iNO4%-o!6a3Zy7ralC()y4KTk(fW@b1jTP*e?PKL$FYrQRg
z5Lme#8M7$lu3-Fd3$23pwnNnUzBz?jLQtN^-k6Tww~c;eNKrd^_RWQ;=#4EKp7uP_
zC696{#>XaNtF=QCFc69088(d^KHwu0?(sTmFSivaQ51Q)s;?mkmuhtN1ow978g*0z
zftbC+f!PEGC9;-n#E#}HsQ+oKv0r38bP*pgDsqmQlaq6NIsVt4fnm?q>a(V9QB^k4
zFaXMgy>8+7BnkC`TVF48
zh#vf`72G5*kR6sZ+Y!izz?p-Yse;;by$&wS=nx3tL&Q^M7yaR1R$$)f!4>e^?eyEO
zLSATjb0QolzPS<7M{|Q8Aq`Amz7?<0pe*HmgAg6RY_)#zJ;KVr(!6c8q*xsb)D~#R
zOe>8D{&@Ioq+{tPwVCecw4qynR_Xz$)r&#khA*XYRaM_EfTJDL*VhHl0u2S_N8Kn*
z;FZ;nCQMExQITDIotbGL3UJ~ky->(RT%)N;6puM#3vKBe
zS;eY6Z0Xz_)@f|IB<&sx@|oM3TnnwYJ=dagEbR!C?Mj4&MdIsKwn!I_Uf^DTR)$&<
zqSOr_&|RIrB7ob9@-itG0dXsY;EGWHFi$m27G`gIb8>J_-+nn8cE0^uDylzz47Ev`
zNs+1br9z-WG4eP6FPar7D+c8UPFMKs_ErZM1b#pc;Y!QEm6;>PxRy@R_r|)F2h+ON
zsUq6TCcm(&A#%JB)6Ydu8-j7NHG6?qteD7QUzaUuy>C>h~6C9fzRm@EMyn
z`{vzG7rfKUDVX(TEvIBbF1@Cyp{J5Kk-F|U%X7kM+@=$`dKubXX4n$>BmvY4_k4+J
z>gPv4&p+WdJ0JNB7vZrxdHqgOv5s+}S)}MS4kDZt)?=58g?RB%nCftyu>22Ek#9UJ
zUm;(u=QBCK83I7w=~~6cOFkH{t6&b9K~q?+Xsp7-YyieL4~*HxcnCa&KcUz;$1
z=)GMb%t`C$>>Jn6#VIJsF*Y|o(oEBN#lE}SwM7PB9IhMzND%2CbS
zB^mhKNhEsM`WLY^xIaY}c7^~*5NaLzUU*^qZY1YMjcZv*mU|zx{D$6xL{h|Cu!}U-
zLIk@yWdx-&lj+vU(m_KfQOk=!4$!M8HbDf4@;~Z0bO2}|Eqe=YoNcif%m&2!jIRX~
z6qLU~>%GtJ=L8Qab5T+1Rkrc2b!~P)&%`ZeaG9C}y|F`0MV8*629Ul=L^xA0L!kab
zRcR{<3c*3P3WK;IJ$BVjyaJ}DKLuM@C8;Te$*%312kS4URs<+=U*_r$-CkZ;j0|4B
z9e|ci<9G4EnP6_6;5NDZZmj*rAMamhZqKM0+&om-#3h>Bf;iJZko_b-Qol>V
z6`7(dAIXb{56n?!K3flJ!VJz4CI&u37^x7G<+g)-i)-BC8WAWe9SAk4fp`W-CC#MD
z3SmEio%1sZw`)(?$?|y#f-U(-N=l{X(Z)S1g4JHJY#CAf9BSlwPWG^$o(lFmduYCI9~0NX%+_pl6c00OVph$#hks{WzAMk58w-46wtqbpTRb>%ps=feS>Pc7v_9r
zk9h}|KrU?*NWOcy7O==Ze%xm+lWn%Oy*gPklWxkk|fdy}--lYwLE8=Nyq&i+r-Xn9z&y@*?T`e$H%BD|p_QIL2(VJK*!xhd2
zDU2-vpNHPkLg38k#a&+X34O}6W9@q!;?V-xiDeozt9l*c^T!JAvGuKG32yW-Z+4t@
zXM`_5WrBo$B*P97k`6|iuK%R&1sPuhh1)Xrrg|YrrPuL1stCaJh@+n_9^&A@4c*Va
zszaM(m}K4K&VH3eSw5aVFQS8swIz?Rp0~zTe2`glJ*$QlPLUEAd%MG
zadmLih%sD@@~Jge3Ze`|42zEwiuIgC$r7kqA52(-3+UjBxR}BwEx6eqxVMD_-rZh`
zucZ;+`2`-mK43(Qwi74dl0*?g9qAy@B5V%#k^
z_iNPNGV)hQ@@#~i&WBX%y`}@MbsIefs=^oo?+S={QY3^icqYgY=(Kf@2^aP(TGFb{
zH1_edyJCk(0m^PHauBX>Gkm7`Uhp=2$*iGrl;2Y=M0y3QIfKnwQ(W
z931uT0^;|qK81_PG?v|j;{(JWR&l(aZnu|f9${ESH@@j2ek$t=yqhT4nzbqrtvxLH
zlt2F0!PbFDGBE3x`u3kAhJ;)MLR2&ho^VyQ}hO%4CydBCPO#KAyFsD9Zo>R20
ziWPo%RVOA_CPvwckdS%$9SQbsCh5=<1;^v@vgzA>AYA6l
zL3&z36s>&PEn-pPXk8;W9DLZZ&iKK_Sa1T$0k-POk`O_?$vl3C89@V#3l9BlO?r24
zbA3n!ICZFc)PthQ#kHkbu!>Z7(h`)mcsg6SRKhjILtVZH(oFMP?x9py1E_mMRO+H@
zJZ=4fYlJAYAvRx<=z{r|Bmy458w-L(V;$ytm`|exeYt`N{l2vnuL(D?&oLe|ZeaGe
z6+AyMfZl#C+CipZC`wM>SS?lHl?~!IBuy;fbuq>7bhZT
z5fK;5jN+E~h$C6!uPnNb85}^e?rOU;b;7Oh&U;!kL~aO=i7A}h8lpi>^EoyjwpR)s|roLVMW_-1Mxk|2RSo7nGDq%3)mG(4kFlve>U
z=XQz7GjkVGW~FYDtfzfdJR^PnOw}cx4ca@Cq8aX_d{jX0q
zAg|iM_suUb{*+@X4lzwI2EQ>)4e!*}sD=`$K#mzmV7*!wlL!qT!OWhwk|RrVK2HT0
z=W~I4pR4Dc`KxilT%pjGJgT!n8UANY>Y{S#Q=wUh$zq+$
zo=^O@mncx`#M~#3k;;I$;wzyM2JTaPW%-_tmL}QV+X<+FtK)l<#?53d3EmL3C?q#K
z%ySP@5B8tsD9A`RUOG`dVoyh?*lyCGh&RkNi%ZIm;S0}wq+^Le+8sTL{;W%hcK7QZ
z6b@RbT#M_Vw+4crId#s5KG7a!DAlLb+qz|{KaWFrq;(*fSXru{&+MbygImZZ6MwB4w*O_-4ap-*)*+LLuyYm~PbBK1}_J9`_
z2V?wbi=ARCJ3JZvgSl_^N`+de<3s@rg^$0FV>91;25NA)<*LjIw|r;GXp=+jBDfC7
zzx-bJZ2%&{{8}=!Q&(~(Y^~?hXu5+X2IRYD%Rxy{C;0@&%8St2buuTL_h{M$i;U>jz01*WYAnQ+
zZYHuW3T~C`QWamlF_3I#=i#Ov?M$uacV0O7cz}t|#h9DRaZ~>gf
zRrsb|CjfdKSg;IzTyCe`*7$3&LJ{D1>eL^Lr#Ba*xJ2>t%XTW@c?9%)tskp|+shry61U?dBu#DfFG}h#Izb8+^5v&_G3c$Z
zbih3rY56XH{bdDJRFZDV0g{<c)`UjFH<}`4~4B^HV4(uU~eOB
zxRJ0R&pFJ@OD_dHud59KBJew!Oig_(gnK4)+wsrwva>cAY|uhChhWPQOdUQWMve&7
zpQXd6C%do4e#~Vmo;hJ6Qqu3JcTs*qqWUR6FXDWL#=hhi&R^9M<%SlU$E&^dYRz9y
z;rf&61-LJe){<+kl1Yy2$`gibeP8~Y)%|6&n${T32=%4yckQuX(jWE+T22dP0?ywi
zp6)Ob+k7FA5O`b>34FNlcp7b2oSoqpdD2~hf{z!xb5K8ke2#_=gaW?f;H~}OOU9Rh
zSEqU@9g4`-E)vCnh5vzyFdBb0I+SdlWxm<%h2`wa4=+P3RlhDCd0tOEj(ou%0p8g%
zq`S^kF6IC1({bgXBbpLis0cfM-Z#v(ddajtYm|_>_^3Tx*B@TyE7!nluHSw7L_~`e
zt<}TC5n-qHSrCM603VGJ0H<#8iT&cf+NhiJd7Hse^#B|9L?>}^WQXq$H%&y6coQKR
zpVkxVIz3%|HLT}wCnpaAK&BX~$3d=(HF%q^$*PPGzTl9>>R;&i9s73HeBDUF{fLfE
zrUDP#)=@8Mps7H>OPXYPI{fIDPuWhyRzR17BCHkwHSmjUPbjghB
zYcHXEYf6GMU-$iuDB3K(7@5SsU)&k{{H4su|JQwBm>4tSXXbFCgKO;m*xXsau0_8$
z`3MHpLt6**a(+*AExr>gw`B~&hOu?;C)=CtZg1X?%uf3r5lek7mK@nGnR~Q;699s}K=)loC@I&sX#z%}kh_>V>JvAoa
zvs^X=WE4KU>LivG>Y7!5TMv0V=W*5>mu*g7)4FTIg6q?HdIkr{x_B5&z2ZhJoqf)uj8;j)QjkeC|>Wpg(
zmn^^0ev0wui?`lMhC2S0kL6#k#oP+=t-gFre^ZE}M`RKzH+0j9bhetbv8AE*Uc>Zr
zhAG0=k{yny?y-|sAUj~iS%3^XvU}BQjFoj%?MFRikPxj<$Uo0+
z$9A7F&NeqkTt);A7mTUlGgWM?`?@ir#P~qg>fk}N=!D7KAz~Kzv)AOO^|RIc5=#!>
z?o@2t{u}XTBZf39_Z`=Us==f|!WTNk-a8$MHluG}7tC}!*MyF->2Y}{g5uN$e6H(s
zKMH;Mv5BDNuC=D5VklYE2?>!HXF9(wV$cZDKppf*l-H9qh{qfcEs;rVB8t139BZ0w
zGO3=9sE-u@CvJ_OFCq+L(f
z+5qDXxXAv1rBv%>g(VAIeZwXi`~9(V%f>fAQ1|K|l3SN4Or?vz#O{Y9OE@~I2F$RTQqKcoWy?ICb{h-D>fs&1T
zp<;JHS~t-b8qa&dtq3_Q7f7r(Z=$i+7OS)tXkKyg8(|zzaW(u2dF{U8L0e*^6H<>;
z++@;5qb8el#tRPh5av_UD}Xa54U6ObW}khGhWVz?$Uo?Rk@b~PResyoG>1NcS?9@k#1>8X`~zJ4r!3?h9lkGUGKxazx%&;jF+!)7_j$>x#pT{?c;gMjtA-h
zb@+5i+2yMnlE6e+nriJKRmH5f8!1U$lC#C*_ZyD(YY2B9%M$m(XFd;Du6qo+O>L+&
z2p>Y}@r)+l6wV%AD>4Wuw|`d@4+&~CGI`VY&|FHqvkSW{0}fu#>wLW8`u
z!LdyD?LkkyovjUhwuvLtmPWJIqu9n;()~ubR^;}KF>>SuCs^K~z;SlBS@J$k(=f>j
z?P3>!m(eP8>qaW=_!sH+7s8zTZ;A?%QqnfgcaJh|OBNabCUdwg4A+Oq3vzQTpwZpdmOg
z=re87J2!`uOl=t>-?(y;(5sn^-`h~>Z@FqB?*S12N`KBek@+BbE9Z~`OHL3~8Vfa7
zMTPIvL?XDQr4MRbO)Ug68YSPzmoH@Tge4#V=j7~o-IrwAg=<;!rL%<`9c^iej#m5m
z=E_pKS=vNAzn9J`M*FgdXW}+rM}~*WazQt|YTpdfgK7V_wRN3lvw=EUC7Gle4L&tC
zgCQ<{08J{#QS<7~buoR}+FB570us<V<8$Tc1lmMN()s-)js{eC%%I
zlhEuipWR1oy>^PsOWAmlKGRWso&02!ZYe7ZUdRog-1}B7p
zb`$wGs^}!9S2)-H?+?0kuCu~R!=s67xIrwKoP?7NB3V`VsKE*aa|^obA78J`S*P^3
zSd9H)Oy1S*9wWX;Y_q;Ud)|pkIY%{J7dicts!~7-Y}p-
z6~1PMzBF3gB)l#jxiUEuJbch5N}1Q&dB&`)RZv_I^^q}{8|gLemroL!>Xr5y=VK@t
z>vV#){I5>CP3ii!^{a~LHtvzp9@flK(jLF6>19QuX6n33)jrWKk1?n{a37vpC|W5)QV&rRiJrkQSV{I(#@FH0lw!!v$`?&JUFt
z{3Z%7D?z1%z%ablOq-!%?R|0Y^mrgkQC@a2Nu+?_b)g1FlcCby=;w}-FBdwuP|2Q2C#AkHxK^vF%gDrhb8Gov0hpcbP4
z`sFNp{`V}Ic-n%F|30gy3qE`)Tk@tM3ncF>R3Mx#3RbTxU6vRaomdnWC7cdX3ztXxO*+`yh8FCaI1xkcW5w
zhNcc{dfcIk4Bl!q8V!9@sY~=4cea%YJ%TtPr@ty1Y%(fdRZ6>D*|a-_r46gj
z%B?@-9hTc?cFGfUNVe5GU4GFguHP}oZdh9sPq1IM&`q`#m^Qk=tM?kw7VK9iMQUAR
zQAf6tJ|fu$(iHpwyq*opY|-I5AJcgliGEnOhUlZk%E0;AOlRuUW=O?g@ZSiRXu~)D
z{d0ru#7jB~{|e@x#B1Mdxq!CVe0RokwjSKtKlPWvIr1p0-n5WvE?y(NY#v+Sp?=ka
z;mv3SPt|(Tn@Pas`L2=hrJ^waEqIX=u2922-Qg$!Dq(%YN-JQO}I1ph2wM@(1lB
z`@$*`t624!VUU&O%N1ISc3xAmvU>Ndw};}D{}$YMo6Lc}_~6mUR7X?u
zzIj_HNoP3BD-Rqx>0>gLI5g0`}b&>ikTkZ!BVlRs@63b9qPP7u;=xb`k#cH%{tK*`t1DZadR8HydCNBge4Ao
zf5%gq?*gxX_-Mv~2*)A{l9A2|&N~ydFjrWUE04KE_6wEDHB$F+{c@Bse8TiI3t0)~
zX4t`BoH9J5wjev&rl1@WP@pK-(bzcYMDW|~?ym+pZ&i9EcZZLMe0R#bu8}2TWxkc8
z9K^UFPuVYBJM3P8pi;!$cIRz0(GR@UD>>i;ZA~gNLdagQk%V1Z^?Ay#A5(Fn-87G~
zU;8_X@XpnvF!aAdvsZ}1xh>xwOu<>5c!5-A&@#&jVYzuTIhkh=Rs4_VuR=eQVU8k^
z@~4^1F7O{Wz?B_9MT0=NLEI!12pSg<8;*0CE0HS5`hA^9TPeoTZ>;KOblo~iegR|K
z5e+n2h36i?vEL*@5Nsq5;EjGb_}JJuGiRy{6TmSXo-*PK5GctGh2lEnbevQODF$~r
z$Y5`I)~de%#20hzMUy0vFRJR>nBtAy?s0Q{7Gh4=`cP0+vn{Z@uhW%&$YT<2?0m|+
zy`_L#1Y!^p*kVbB8Owgw{qgKO{6TJa)uRSF)R5*d0uBKkE^c2%}Y?YqbLCdtndE($ScrM!s~nGx-{GpQ{2gydw39KwlP
zZxi1tYtF2c=OgH#38LKAydl$4ZFbbJP&-f3IZM_95$=|7{Kd|ew??pp?vr3qiMp>8
ztnIqDJF1^e^>bc+h*s34@Vl)FX7e9!n^eB?PhyzEkpI%}(vV|(6rWMy#w$>k2pG2#;m`$qeB`quhw-7D7oCl3)Q^ziz#+q?!v10ziE`j0g@n1AdtQcl+89;S%$
zn`o7-5nB128%Jv*9A)?(p-prrW(P+Pr|iJb5d7oUm&LMD^rfI2cTy)yX@nnw1@mjWJ#xPo?(b^I{M-s+GKlj;`VFhOgxC
zUQ$vcJO=`4=GK`CZ+vk%t}^Dlqypb(;!TNn==E<8srcf7T~or#i=^v<4N*<$-dLW{
zM(iPDyKrievCx}ikhulYVU&}m=hFTDXhBeS;jYJ~mz#uX5-Xh2Nl{{?cw~jgg>j!-
z({iUA(0s_9lwtP~Dca5z{_etRqjol%!fWad`QozUtO}#ca*^{H7Y+%@$P&cH8Q{9O
zztSrbl$V!ZSgyu8IeFFM`jOJu5{k051zubn-h+SfXnRh
z$Z;WGGP;7e_mDr1pF%&EE*vaNb4q%1pi>f?KK*dP@GnYM{iqPZ@wYx*yo`oCqrl``
zw0J5$vT(ZCvJ^#~%#?F;k6i^sEs|YaoVecWsgj=m7NUeLDcen~of&bsIEollle-EF
zwgQtM5SE?}lU%Bkk~vMJF@I72?c@gj
zRz5Vb-gTJR=9{5xGxPo&7Q?lwA%i;mIEA+0tTLwlya)X3Hse<5>et!hsL|1~^E>`=
z8c7hF;X6ir*8{e`P1egM{7as{z8?z4<7K~Oyf62cyttyaZ9-Xhvwz_IsHiJOA(iwQ
zpMW4!&-Ta^9q4Lg#VtCUbpuc*QZ5_~=HEp!;J*xXBc$fuZVB=<0^>SawVy+#KDYR@
z5mV9bQc5t!n>)B+QqlHUqMy#S5GeX~(-}3^^zDA3+p*sNRA|xSW%$8~6b3AgBlUIf
zlDNmRX2yBIfK+b`sp{74J=eaaffg@Q-_|^MC@#UY&IA3G@{(s
z(Bm#G#lZeQc)eTjXS*R!k;61E&5KmC^9@+;W3O=WCg0~Qrh`M7p8eOb(HkbZ8TW5-GB4foCRhlYlo^P+WhyF@Fd;alw}F8SkRL3yl-A`yH}pr)p=>!t5=4dYGE^&!t2
zTPf~jYyyyCEY!gG?G50h7?0`)|EXFWWY~I);>acwWj++Mq9E=bx@?KLI#C}VW$(>Uc%>?ZCKpKy6ihfTLp@OZ9r;Jv7)rcCZHZ>s()&zI&
zyHn@dH4I!;5=RLp`)Nzq=}Imb&O1`Fs!4Qwlu=T3{_HoA8Ta;7v4HdE-7i9fU37%j
z4nn`M15K)*)#&O0TuC8VmPU(ZU~x???rcvvppXsiu5Q8`K;330iD~
zFPV2qJ?KN)cs?iS&X-TJ^Mq$7>b{)St<--f!T}cbgsD(?D9wS1EooJ$@Dk`eXZ5%?
z@jk=dh$M5S_GQz!*gmh3xy1&Vjw1%433pv1@f1@xubL%eC4)Q)2;Ah$w-~i|I-_iKiYz5@
zQI7f?fG>bS*}t|wH>B824yz)YdlI7dq1J#3Keyik8WNXHt^rUwF<~LucfX2Y`7tR&
zT>E^*q!?s3G55GKo?zvl`$Gy8^R<4tl^p(xJ>_ae4kgI~`#aLd;E?26mc-*7qJ;;w
zBaglUp^JsR7ipaV2WFCxdulw8=x!T9Ve!lmD5=6JaV^h?Bo=tANy;4cmkqcPf~^M4
z7}`z=I^SwWlj}hcpWoJNX+B+DnZFh#Tp)=yTl&`CsvvBah6z&#*M1Hr;ih~q&6bRQ
zcQBuyY(&V7{-A+?|6^N5SJaMWt%u~$OW-RiQ2!-mJmE3ocQYBlHj`(Kbc*kM*tI48J~6b;9#-I~)5p#3GLyNQGMvWq+Kr
zyuEAc5u;KCxi2o7W<>BtMUYfBtF_V_ICc|GzW|GK`@>G_B+mzaLGTXK(b}BbpfZ1z
z-=FUYV!-ZeaBT0Ey6)G5f$8X6gZVP)6)uvtUerQfjrK!HCILc@_9@DwClV$gzg6B7
z7VNJCg+cVCvyt+Piyltlbv3THsll}*7Dkp04LWZNuU=H?CWj#hS$#ujjHALOG*1+o
zu}jp?7Z?NYpb@i|#v_ZB`X
zQQ6;Kr9qAR@4B~!`yl1D;$BuLzoz4b;Y&wr4Q~0vC=(^n(dTfUPaPx`u`Vt+>qD<`
z8%gN{clJJOFz}K_eh(}Dzp6yfk-L0ROt&@dQJVadw-wSTn(Jj5&&i?3lgQm~>5~BI
zH)%O~+W9SGQg7mO)C{T);e7E7m^apQYoqbb^>ZI4NlAThZ18j-NFKA)LKHa|16rp{
zd9xQ-A*0evQSMyB{WShwiXMcY=*{u;8$)y~|F{i6++pUV_3Fwo20qj%yocAx4WRr4
z%y1kbf{b-xSm78j7%>qz;ZP?i^YlB;MNy+qLuHMj==@MT{lT!?A7s`Nv?&@2~Z|8M9Njw1`1*BG&(XfdQCazccXt^Uk{_eh=amzo$G8G
zek9*k4>5nN)Jmy;f|1273c1ZA3h<2-=SiPE8#RE{kxSSz`Yg*jp;{o66J2@j8;WuJ
zOa136&wEZMN0axKYuB*37J)weWcRbLF@477dSDpw0%AhH!OkFOlB4YDfO9i{%j=#-
zl_H!7o(L5W{tiotorHUfdgGmNO~Y&%n+i{ieFl)rDzH8ZT^VU
zjP!v;(c&!mBmD5aEeSBw8&l10e6P)@e2(tF23p+>>
z%X^`y8W8%nm%9Jsbqx(5swBe{Dtt4-JzZ%$?VMJpXsQewnzX7l?#%*}2erxO(IYat
z2BE5CI0zCqJ3FiFN?x*llDVxz&@ON|LyF49vsAV;6&l_$j(QvdHXrU&T5aVSE2a^3
zbTyxPYMSR@#%dxbIp7q1eY+-?@t;(Mq>|g
zJICsfBxAg$_Cshb_g~6VU>2FiH_7^Ey9|;D>0%6B;@VY8mpc$ur+S9w
zTt>sWAs>Ft?lEh
zHn3!0`qKd6(ZrLq&?Sa0N$Zmu34m30`*)$q@DMvX5=U_
zb=Z^1IJy#^5DStmm0T(4`-FCr-X{LQW0p6zfrcMR!9*P)jVJD*5hJ%7zxni*wp`x&
zsGwF&FGLcJrZL)0I2pNX~CJ
zOq(}cREb7F9N65tHxAL(smk}#(S>dE5}-yubR(Q(UR8ucBO`PE>2=SlC=|%&j9zCe
zwQ)Yv%*w_vw5fO7*5wfTKPX%-!07RW5G2MQpQUDCV*W&mYtfUOF4Xm~hOWzPpzto1
z4qwfL%I^L?_7RsT^=xt{{7`%J6Kce%dub`W45YQ)9R}}vP+&@L@|U_e@0$AejI6E;
znQ^74Kx>AzF9;=5j6p2-a88&ZQ==9aX2yADcZe~LpmN;s;Ne(XEXuDM4ORyQ={uE;
zw!F|p1Z^|^bCBAA`5u=BkK!;|0^^ULE_c>san7YZZ>{?mM9x5WWxo&TTROxxO6xmU
zGKZf=hd8BNEi5MOS+H+pG(^Hu{CqWXL`rSyC{p~sb+fXjGgHEF2z=-B{+sfEfd#!n
zTU!Zb3#o-1BKTHYtoC680R7{;V6EyNR38THd20@OAcDMO{BM+Xl_5MJGcMvBt&9!k
z|2Qq=-}6ZkF;8!f(amWD=f5KBL9>w&o5NPNNWF!^)qGSMlq*
zMST!1Zxb%Iph+o4QBzF$myB!{Ejpq)jc9M}J*XtKt>fCdkBZ%&^|W)6bwb&+^0SnE
zi~0S%O>~C_j4Kj*YyR8fLTR`>F1jt%iwzni0f`B}9xTc^3+eYb?XmnHV7r8Xh|hM2dTJi{IjTrmRgdKDFRHe9nOcNM|t4U
z{y%@hfb`tlnu-`7XGtd+rgmI=ahF2HG^Htwwbp-Md#2H&VvGMkop7jH==QptkzlF5LEy@Jt`!LO(CHOE$}e5>fOrv(PMGzUyW(4CVWRrKNZ=V9)p|
zsJDkL-v1+#!|i4gtK_m54?uBue@bk8@(EtY`5&GjF`&*fMju+9rKX+zRx`I{E`lOz
z|7UW>+ETDj$P965=?)hmHGY=7ToB(=RLC<#?&)XMlx4&WF`xiad}V*;iWQ8BZ#if=
zgTeICj6N4XFKpw&B+QZCN|A~7&Cwu_rK#K#%`{#MlEeVVjc~w1%2TrB4T-K5CF|)1
z%m*3yVfvX*qpT?@8(F)`3JFlY7D@&N`THMKC-QpAEGU9eOf&AytFjqn68^sG^VhkW
zhLD|IvF0tVkjs5zZ*ri2L`D`e+e%Q@*psjU|H4ZHhn7^|f)O*wHLE<}Ek6@iwL53X
zXsn?~$cmIxivjmx9^^tu1rF^ns(PDhI0eCIjhbu1|8m6)Olts;TZ;weg+6_U^1sdo
z7N3A}iL$_t-VDW(|M5i5s#=GwogqXw##>C#`aWyuTPzWuPuB?yP0R2HZ8j^)Q7
zL7xRwD{1)4i4nOXXr@#Z!IIz=`eHIz{kjO!Z$y=wqb`&g
zoo$aamhnKM4YP>v)`3nu*+5HNg}oM1_HCBJ8?JO92!(j5LJ)Xbao^)|yffpu3;2-5
zn?SI}4b)C!pQmhfYMFhv!r`#EwSQY=dxjulM1dre*A)!3+&Uv*0B}e2KPNFCgE|q8
zZfmqQ#VZrl)~^KNclCkP>N4%=0xjlJT=4_t(u8AA4kr?;a
zf0{0X^hQ`+hzSxA{KnAn_hAqmu-fQw*_P4aKe?7rZ`Nmr_n
zN8Nm^J8VtlWz*PWYza*}ps$;&v6@7rc;1wB6f%6in{`_YrAt56W*WUVxNTL;Qp!}K
zX9WZaJrN{Cq@>-~kNDT4tFG_-DgP)f>s*OE@}^%9=ldi(VcMXlUwlBe8koj0wjwXZ
zPe1v9Oaz%Kp#FTxxG5Xjj#W547#^jW%o~$k-8D9ETR3FZS=CWs^_MR2>%(=f7-;_t
z%U=gU`2dYRfNZ&d>n43s+Nz`bW0(>ekmbJ}VvFg}Aks_x6^7B)&e0vA4h#k>4KW}E
zSC!1q$!kiNC5z2&_xsbCv`sgZtEbjCU2i)wVcZ8)Y+lmP05y8bIt(t2=qSYON%?bL
z63;e@sB8M+svQxArOr2)QXh!!5>F+CsOfK7`_=oaCB{$RFm-;*z+m>iZvS_i>=YtK
zO@7?8|8D7Srd~y{N{vlCq-TBpYIfV9JfU<;V|%{0=BcYL5DS=U@mxhSyomtU`rOz1
z48jBB%ACilj)%a{YmT0$pXVy-neERuWO{mg-ELNc
z`5e|HCFX3GTOTz&Rug?U+}qC6j+Kh$i1sHZCxKkFys4UThbJc|2ZJEiF0gcpDSR5L
z$jZs-v?Lk^_SY%xY5cN1*XARFC3FE6x$RfvwO@f*e!8l^UB&Y0Ll!=P6$VuN)#k$~
zhQuWfyzNg*&Un9xjtqN?s_Ax)iKxTGbK_i%$rxVyWprM{SM8EcoP$C6<_a7fJW{#T
z1VE<+(W*(ny`P1dy&jKUA0?R=iexOA<&;1!J0vZo5Z{zkY_>BLfHxWDDGMd=ZEzNf
z4$?Ck6Ia4QRDDg6hdB8uS_S;s?>*~Z@`-VCr3F^CsSVZWt1O-1l7dEK{rH$N6q$1v_FIx7Ci;D-wcf+R;1z}AiXg*WJ^LPepoiXyR8q%`js?5*ZS0F(VdH%i?
zH|~)7hGyS5^6zqAruMns&)OPi+66o`((Q2d(pwX*nAr&&iuGNc5d|egRVxkKKWMdg
z#deVn*-{(!ReZ_BxxVJ9wOb+DDkd*@1s}O&pqWjIKu2jrn&wqEaYR4c}uLpqXk2jVN8~yH}b1jY*8YZc{nDcj_XeQHsRsuGT&>vah5g|o&ca2&#zU0or?WmVb2Y-EC7#F6BI
zPu_joZ0qYG0C=^)T00xe(sELLVa5ag&5W;qQlyS+dpxw`U}NhD!XUZX_Ia`}I8T|S
z3~(J`_x4pSP;fX|5esK?O_%oY5YW_fKvWsn#;d*cOxJ9?n|_!!e8hj*ShKsoZ^}zv
z^p#3nzqf2*bh_-s`d&}=hl-D{X>TO(>eamKA{089!UcKPuxM;;9o}~Hz0JtP#J}9<
zabHroo6O~}*t^7a!{QuQ|onwvX;F!Rr7-8CYXHeU3dFUMK;l|gl&s8~*wq}0M|zjP4#>)8y^W91wnvg<_u1dB=RcP{L(AGvfBJ*{N3O`r
z7ThqfguE%%af9a?HyKfHspDmnH_v`N^2u40>bCh6C0D2Nd-(e;8JGF|A<+`3>3Cwc
zTb7DrnE0}}=}H5^MnfaS&UBP=((lSNQS?Z@Tt(?2fIJj(4lr=FZy89vV;xFFR8+Lk
zE|duf*hjLMpROLVn(v>iW*A98AeP`eOVZ+0qwdfn`W0hTqkA4y`n3-IzR!B_Sef63
zqOD}PMSYAe=e#E9$hLdm02vu;A^etQaHgJk&&mcKQ+%#mCi!`9J+)wZyEz{2wi!d1
znYn!}d$C@nQt0ueDuoOX$$BK+a-P8mY5&}eI$r;iY*CN$iz^WNNynh5>
zEFR%(q(bg-_s4A+B(@EH%Xg*AaOrI~$k1&}4JFG`A`%kQ#ip;Lx=N5z6n|gfPP6H$
z69F<3b`x<
zzniJNxoXwRo&g?K+xdGyluo|f8?%nQd$hB;
zu*g$wUU4gM6cMo+5w69iQ*4Y$>PsViTGOFQynY)}NA1;flJFoXagOW;JS!1AT59q4
z*f@*WzKnQmch(h}qS{=0OuH|DFPFxZ&YQT}XOgh{OS+L6R*VwERO(O9NXX#b9{OcL
z|HhZAY@EvU2MmE_kdW5x?XO=jSX`7B&Uns)1(k9sDmpUl-HH!y*Q
z7>xsKrUE;AA$(F-+ZSeTIf!pvL&nSFP5`@n)H4x?VV7j%vO%ji@>y!dS*S12gE><%)H@80?n~ur1w{CZLcLq}nkWJOKpk^|?|FxM-b}5N%g8>pfKPV0Zl;(KP>rKX-#!N0Y!YK4lsl&(s^|4pHI)Ru
zr+UN{4J|3}HIy!6kSpy-6YZ;MN;*t+b<&sZOdTn7bgi5=v&ET1-DQYAYT{FR8)4j5QeL3ayy2iC&Ef@>GaQv1u^PW#2Lv*Yvg
zjh*(#^B?k5{yI&{|8-rD@LYkXco_jlHN0P#o%`*=a_sYbe>>MJ-19P5_FV@hGj7cf
znH(mi6?xXs=Mhm;ujWMu<($wJj8)#S)!jV4+5l$XWhR-8th61bSH|P^6v~aeNncJv
zIFak$#;*TsVln+Z`fL@kw*Y19r&ludO;~1Hjln{9Q~nm`o2F}YNI5hLJ-d_$2GYJN
zGj(~>-b3W$n~{g!c}H9X4UmUH4qTQpHOuQ%{=nPYQ5%NGgA)Dq{GVW+MTEw`1N+)G
zBiGin+bX3dn53cr6A6H9>#bH5E8AKAG_y;>_6;{CovsO-e1m#?0*@=o&s@A;Q6ot!
z$&q#I%Sy~b+Uz+$wUiC_^tE=07n;!e9kRdEzq-J7fw$Ti=D&=!Kqpe?bCUIGEqXUwT=S2?1=C$|n
z_0G|efkVkQQ(+)DJTmfZJjy$+zU9!y^Qt;y6JmHz{IYQqxRvH(uVW8`eL)Iy3ya44
zqu+j8Ty3|*b_U9P4pg*|j~_pB_OQHBCptVlytuk*I`ip8qzP!h_VhePeB#iMB_yrs
z)P*!WICd}8YWU`>%-|{R$BUfn(EPCWdSlA5;N?}(Y%M*9flg`LjLpV-^AY*E`{5l5
z$kK*m(_&QTTeEc?!>|3XFP!YLXYW;Kxo)uQCX-Fk!%pFF934de!HT1TrUhY_vvh*17nJO|bWFFltdT_0g0iSQ%J~J@zJE_rEdcM&N^t+KdJwU)
zOG=frv?*u%m*%^n`yT*jP)9r8;&Id%tl6@s?&RWPctTvDkhV5eq#6;<-CjFu0Cn}E
zibuf3#dS^Z-5UOxQE9k#_WXy%J4WER88;8_WfNCFkQa_#r5@S}1s&fX{B4Y70qBNY
z6^7P_O_p;22S!Jb@HnnV@EF8d-<_@l6#Evsifqk{S@ry@>F12p=S7$0#fDiG**FzF
z;EJ7`oSJvO0B_cS&(D*j$^AwY0D**jEdMBke!Y(Ie}wC{05-I4KD_N1TC3|bs#}bd
zIj>;fdIe)CJZg?j!m|(-Ap6H#{a0lhLt>VOgk;R{UilMAG6G0P|FLVV4%W9tvviis
zgDA;tId?z|GE`}rL?aFE>f(k0bY$PCM@Xku8_fRiM1fj_FzBO#!dd!tbMVK%86H%$
zTOOe|$Gf|$3y^ZecSACUAhBAVPe4L1(XF@rW6?vH8%N4(@B4C?9p{3==gFfdOV_Q}
z+3`DQR40N_S69wY2cu8nK&|3%z32jb+cStvz`5CyLx>6iOYj7ymKT>kz-Tbxt)AOv
ze8lW4V4Q4o^V3nodL;gFeQmA7erY{-D_t3C`)P3xsjzk6Zn2>aqI1y{UQI3W)`rZ_
zI
zq!&PacnWQcASA2BN0AYwo*%tFyf<1_9I^WPT2$^$o5O)whDs5`gKjSZ<22#2#Sa_J
z^2sLWDqHXzw=X>M!Q+tG%7zTd#$40%(%_sCcAk+5-1)_15{JbB?YBxc5)ct##Uw^5
zmT#9szx^fGzlMpl*WJi7o-PU+fyBuIEC)qb#;f#se!@p`+Mg}Ac)YuOMLQJh;QUf%
z_&m<~$Y&|>9=8He<+*Qa3Xf;e376uwPx@Vxs>s}ZrNjFB7v5z3UfARMXtM*E!gIVg
zG&{uW{t7$DSR%=8xs`1BI%BeYWuCR?q&nj<{2MYBHzQun%)jV_U&Uo`sW&do2ywDW
z$6lieUG??Ql6o*<`_QPV_bJ*?PL#OyoW+j=A8LQo)2KM^Z-&M5*ZW_8FCLWOZxl^&
z$j$5eB-R8u5>;0iVZdd0;t;$l6%m`@$Uwl!oEKlY*Px`^Jex1^I=a{Jr0ejT-tzDF
zkm>iKeu`as2w(*!ZqggySEZ>qN5~=$eEQu2v_<`TDe9jPtE*>V?49rgvHUK@5dbZ7
zQjT<1IaVd)m5KF`Fuy?BjSKHPwtuvuERftdJ?HsPubA;_XNHD{UnanDYQy6esx7cu
zp9RA^M6GsRq*JZ0i-5xw$zpr1><$ZbR{D
zL($-*3A6Y6A5x((KKH4F4mp(gFfN1(!J*wrZb*wMvWe?l%4OapOr4?-3KGnF2*9Y^
z9J$zt8e0#UPiMHJ+BT1eEo+P2n78cuR
zWsJbc6Em<%<$$+k@ptKs0(m7XzL1bc_?uYawG1Cz5k>ZSagbIUxkDGm%cy%QKVoYX
z;J7sRILKbE`AIgzeJm`dzto`~J~ynucDA~l(*WN;QoJshFY3=Ol9V<6Q%1|f=O}7u
zP;}gR`MY_;w7IWA
zEDz{IO95zC84p{}gr*o$)`;Z0yKB{{QRx6nd%bmYx(=`hoWEX}+V^&L$A~1a-$0kY
zj0*jAhT24MT_4PMZ3|!Xl{KFY$mtNR+CQie0K}D`KOl|Ud=B(~D;-5G2{sZOFq{p)!lrTX`(mTIJ
zaKo`?&7aSXMj(%cgho8;^irg#h4my!F<{fsB!!|3P2pwh)gk$o9ufde*&i$a7CGER
zd;%&ez)`s*kNWyRK!w$~1E&V+SX-lpnYT`G;(uhqZm_YjfdP0EC&9V{akd@+L7w(f
zSU5OTj>A%>!e#Tet)W_0Z~*D0QPO&w%v
zaGxU+QLAgr0*@AILzWCaHipkt{i%azW`Q
zrsJfi5{M)k7>*2e5iJc$&Ik?++|sH1;PZah$B`ngs=RBAoGv*a8i6mn&}~s>P`&EC
z=Lj_$+DGC2y{acfnRaSE4aiC-0NtJbf$nXU4-eZSb|AGZBI*cMkuP~3}mPA%`
zB0|D{k^L_u(mOYr+VW?Pj
zpaZ@$^<|h$8q2)w6bg2#v;x>2JXH1=>abOjVtdNK=14ZhL-uFCaScY{*RD
zW6D0plDq7~B@QAWAb2T0$8~7(BxU=7%Kruc*-s_YI6%D#Pyv9VZ+(3oqyG8n$-As^
zi}SmH2!stgFq-V{03LqXE{QR8#60ars9Z2z@f=;>$1N5)
z7u01_rWXazc`)XAv7i9sC6RLe|Io=iKj<<%tHu+{Yv24FuKK#KG3(<^es#16ZP1iu
z=w#ivg4%752poIP4EX1M&1`8&TUs(OI1ST03HG@OwuaV4l3^>x1jDVnS|n1F2D1X|Jiv
z2*J0_!w=V&^L4@E9U8_BY%k{U&%oj-iPS{b)3BX40OG|#H8W0MzdB5-Yx1!zh*KXu
zpNl+KQ+*C+)B8HD4(zZuq>Akq+$L^42C!0*TG+=o>ymBjo8Q6C8-G_a4}A=E#GPDp
zuOBx*F{wZojLZ~FXL
zf4*&-x0fmN((~w2@?*Dpy~D$O>bxit5}{6zjn0A!BO>B%cJ0RWuiE!HY$*+YTqmD9
z&{`(DnOIuRJG(hg#YrM~-i5w5<+SBg3Cjy0%MxXP!c?c7d$6RR3vBb4UQXZCQ3__5
z!KLy+AM#aQ>y9<`+m4Wr65ok=T@fwse;M^7CV4#HbbR>mX4UhNi#`UXwg%5&J&NCP
zD7U=_iB)){^IQEQQtqTuBkA}1-Fwz1B9O0ao~tJ{G2)uTGF5xCf$(Vmj~_UNhQNZd
zu5=;KkDp$b*HeOc@MT~?{*)jnDz&(9MsNMH_x>X~3{X4H0o11X(n5md=n?#rPHk;1
zRr1HDQd^XXw_`nS2*!ih5&q8(EsGIf9F)R-8CwU%0%M_jW;F}C!Y1Oq
zJ#~VLFST$W?-j)$a$oj^yLilyTV%?#t_B{m|4D=Tlksb&@=cE6Z+n5ksLW{>0JnrSVLzSS#h!4f7b0!LObxieig#}D$X)D
zcg+rb-`(WTQ&Pe(RQ${067!))7>-v`tz$Q<|R$t
zZL2}SNpnSt^5vN@Zk|+Ubn9yFoLHDW3$NHU07m(7MV2YZil9@
zh~ZdoYNlrqHlN;zenaK_I
z2(80L4YUWI^RT`uwc(3X+1cOk+S_CF@1V4XzUL+FsWk27zMnx@-NFv_1w<9Vl>lL<
z-TU^F5-nWQxq^~3GSgA-`ro_sVgDa8g_Ue&{=c)`4UEc*ugqEPMyB>*v^3{yGcyUA
znCce4kPcRi9cc5|p{N$mC<{L6LKOYKE
z(H}OKeM5|(b+b?9#m|e2i-pC-xRsK|9=NEmWMmafuJAL{#d5
ztDPkevFkKL2t
zVLzJ%;6Hk}``)=PS$sx~8}^SRg>7o`{l%DB7WwVaM9RCm#c&^nndQsZWd|)&BX4pw
zbwqr0|AYh>p;3y^n>0{-JR*Ug+27yf9w#_(35Y|ZUAxEI+XDlcGc&yNV9am)0%B2u
zo$F{g4}h&p{bTDu&uFzX5Djtf5~!2t-@d)HK0?5}dG+em>3Wn10VpX=&_FM1O{>TD0LiQb3o_Srt%rd4?H@xOK&wN~9adhd0
zS9D)IuWjV`_hMGD$PTY5Kl47oRhvLGR
zU9}|v4ep`QQ4a`ysHB;i_+0YN?R{6UU;Q$YznIqZEmp~P#MGmBi4?GJ!#$2pxF(ter>%?;m4rE=lAb{>{z$Y^foX|61Xi<_n>7b
zQTyY(cH(h|2Gh(s0(zG1uMzL(1G5T(V;R-ec^+->5mvAN@my|8*mj_YavU{2Fdz%G
z!PImX?Ar6>5-R|ux$n26g#{g;J>SQVB0w-bi>&uMsGX@C#%1x3SfkAsaGF9_ayrtc
zV=JyG<;p%x(l%QUoKq)EON&z56Dt4LFfOEz>RzfmMm7-{as_U`A+UvMNP@*3CUmf5
zt}HCo*{up6z^9E>+DD{?MXR%1g6NyJ=Rfdt>Q^EKc-aU+Cm+P$tZUy;;haj7f7FL9
zMa%33n_Fa=TZmm=dPZ8}IbAObE>U*0QCbuKb5V%040Cl69UWo-14o~i8Q0nI4H^6U
zA6i(W%s6rWw4Rj}%8A{Y`j*;S{MWBvvolg1X#GE?t}?31E!xsj(g-3A5`vUSHxkle
zAks*e(hbrfAcAzaN=SD|3koPmcL+*JH@x-TH{KZU2gBZeVR^jh`^P}>?6W&()_n=|5RDISTHLNe}XgSdY;zw7b1eQaICgUnvtVV`nY8C
zZiH(#%Z}&GQ#cw~yz{!WCLTAFQcP^RCUrGr(cLtwsz7)4HVpWGA18W=K3#1IPFSyb
zCCqqJP`GANo`i9Q6}@C+WrcHbaq*J6=lLg~#u!M#dv)O%q19cNptZ0Ow@Xe`+S$8C
zfU>@x^;1j8)Ko84F#4B!j(b9v-KD9ilrG$Nh2SLEd9lZ13wzkXfB)hZALU7`Qqb4f5W_VJ;Sm6gpn^?Udb2jl?_#>=0W
zItjt~`L}^xcK_o0?Agu9sVQQ7tlOyc8+za@^7MMPwmrgLQ=XcN-rcis{tvxabPHv3
zrr&A*q-h&1i{tmp;D7r}MDuN4R^h8#{e437hR^pJaNnwV1q3uSXWHknYMj^eqd95E
z8a}^%XX>yK>)VBkk%K*lMA@Gw8XJs8mN(d(S|f0artFuVU6pca$YDO>>D5OmDarb6
zzEw(uVXPgC;ZKhp@UKjv_K1*>&=hhwJQ7NLM%tL%=H?k?Z0t*2{RU?oQ&Uq7r7snn
z`~m`tZLiT&eY^rbf6#CfEH2ClaJm$;w8f8CY&+Q!EmD^ZPDw?-O_DXBW@Lt0{q9>PAodGT>G|1^7`Ya*P~C0IGzewX+6V<=+eP{dF5NfPKIGT
z{J$UMv)$U}W>;LB>9r@QPglp#>SXp7294>V{|4yAOC=X
zmPKLU@A8j+L|sx^zm>q+1|43(ZEhvuu!2vW&i`-|y`*&a&BY|u)MWCP>KjD1R;b$P;+`M@LV^d%qoKPZ<+bO
zvLitfjEJzAYj+@D=`XoNcKQFc)?oaSO9ZDu1U6BSstm^0%k;G))!EjPuRU~r8$R1Y
z1oo{LuKf4!Ka{e`x|;Fb>7p<=phJPZMHJH`cdkyviRXMSbT-l`v^pRSOFpSfkfd)T
z_5Q2B3GZQDW@hL;NtgyHqUcHoSt-MtWS|$@(i->g
z@(zhW_p0|SUsP)KNNXwZkRbJ$?!%3_@#rx*Ed{>gs;tmAwD8DCTYbu?fWAKW
zW(}G+^SwcobNbr)`kP;jhDaOk(TN5xR%3p*pHwFJ@AZa-02{;2|rkfrqC-N~*E6
z#5qENO*U|OrAcpJwazi~5w_DYOKh0KwKW#%dDkm4y^eG1CnpDodeNdZH*VktUTK7(
zU%gBszg%CRIOR#j1Tc&(#eaG(jC-)v}O!?D$R8EH%}b8{?hZSA7s;s-+1T<)2!5tGe!Q@xE&
z4HLFEDl049H><}!x3p|E`!WVB34yMaVPry2L?kUG{Zg}|BkR!hu*}flca~t{Qd>jC
zySr$2@s#j3$nk#Hk9wli-D$BPWJU2kiwayVU@_V@m>sii<-R<>DM@iI&zUhN*g7^g
z@M4_grE|RGn5wFjhsQ!%s~KF0ora>wjN6FINIciBF=rw4Z8RW1Uy2d+4tJ=o&b+<<
zq;6IxtPHg6(d&ej?!Zd;@ZQ6
zhUt&hCGc@d%BC`BXAL<@YHM4)J~@|{kXx1L+BM%9#hfO_z?^Vvr+W7`c;oAYMaCRm
zTrzXT*$7m0ABU|o$-|QSrSJj)!bQnC|
zOz%jpsRlWqC;yz=r`
z&OHI@=Bi$=8B&^VGzw6B>Cmyk{~EX
zRo#H;-kk97dZ;>$tHV|OOc<{zr-O%f+zO5Q!*Z=_S^2bW;k3u2ReZC;Wp*91HcTk=
z8&5K;s}sIt7hnyqSHH>szIkqI8c2jznS-XJMD&Gzc{z=9cE;Z3Rs%tUy~o*%J;%?i
z4>%t`m3+JsoV6que#EYWE0-yA6XldWlz?WD(C--C#l^+F*}%lq^i|_bA(?YZfiWf4
z7bBw9xtOtaL~?my#QntPz1HH`mEFhq}UkP
zSb_MIJpKqB_vXzT?cbG?Jquu;^m|rrf@x5en3P0mVq!9+`*f*mH{-s50E`K&I~aYQ
zEbL?+df3H-PXMVI1T}H65O%}clJ%CDdq!7u3u4>;uzEV9zJSG^KqeRGg|kO5}7fG8DPeGxuRor(Ck5&48Y$8&q;>%93j+-m?KR4nf!8K`&to<&A(dg=3p$Yzl5HHEp&|d`G6nk?
z>uP*X4%unf)D=y_V9p=62x_nA*UBD{l9nXCbkvkMb1J;&E_jP6uawEeL|0-?WM)2S
zP@KS+QF2i`7a5^xh~E!PHU}02RBSZjisL5
z43e>Ys%I>RSCqmCJ3GP5o8ag1TK9gIk{lczc$+);U@PbBxu1UxX}@WLyA1whwAE*xs={p^{dB_
zD{z8kAa9!r9PI|*Som|`x#QV<>@88Cn%e#5=Do>_mkuc`2L}h|$32n-aB|K16RLM)
zGH|^bK~1GYW}Q#E7k(kBb7k8L-+QY?f{&$6(k&CFO(kf{KokBz&2xxC=52cV_kQly
zZ}9s(?^3nR9)a5BF9^k}za@Xf+;-3QQRe>1^78V%L+?Z_k|=q2NAMsSf{fbyQU9}h
zc&?YIr0|E2kB{Fo5KJkEefO>%jK6s1gBjn!p^xr+HZ}L|?b|G|eDRp8tjwg5$xe9U!_FQSPFmtIF@vL0X@~1?Opira6G6_xOdlCZd~lGk!z9tt
zwKqq78WZ1ZzQcWoBNm1vvZzSuxx6R2!j+?w>z$2;$p+X!iGdi7GR%(?+VSY=(@N(*
zS_#~STVg{G*^`A{F@i;a!nFGh0|Nsp`)vRdb`jDCHx^zd0j)-YjS(iFJ2E=DyDE8=
z433!@;!j#iVBPPzIGP)TNDSCS7F*iWQ$$!;*kUFW(rI6RKRXA@@rq8BZklVrd3M!Er+Hr5#$_co9<$?T&?GhT#Xxii;Ihs
zX@2atRaI3VJ$|gIsrh9fU33D3Q~)-K{c;8dwA`)dd?zO-&37>J;^Pw%Aa0=Pbj?bf
zZQS!c_++U0`0`k#aIPe22s(os(>LL{J-n#?Ck7_%N+?iU=a|yyMzG5t;TH$
zcu#|CU#olkz@l4!b$RBomK}OLTH6xbv!wZWgj+*p_{CVMJ`!Xhz&e@ZPlEXmjwCF=
z9rSpaAF{pQQ)9PV@bdF}cO5^W7#Yx-jIb@gMo)^Q$OVTHeOz40nZ|Hd)AVVvzW?B1
z{53QU!#S<6vg&zBeSUu{-GxIlvtAm};uO^Ia1Bl4m0ecW=X7q3&mH%fi2KhRQl3!y
z@lR@zwh5p}DfhkJToA-xo~(WwAJp3?Ln{Fe8m5kq|1#rWJU*m&qAYbYT*@@}ATcr2
z6$bbDzB9X$-H!Y*U<=rO%i3
zr5$%lQCHhkvulL3%LuCC!B-h(W?*;?6+fTL-w^?i)|$d~+c%WF0Sg_WOhHDM=O>gQ
zU*LM;sGg+-xDLT&w#kE8@G2>_z~jkwa+G)shyuHCW6WN4nNJjm-l9&TQ=V|4DAliD
zzo1i{&kPrA%93z4O(7N>d$$Y^bf~Smuco`D^}YtX^Dg4XgcKqiP;*dUBkl`(xuSv(
zqF*1Wx|lRS)ork+<3w$%1`**?3%xBa=Ct5SPEKyx|Ec`w$rH5EgI%}{!wKv)5OfCJ
z(>85Pf8@xKQBa7=r$7*eTUqA*bx{~1x7wuf;&!*Nl(t@a&ibC&C2&4#PZM@JRaL9c
zS+Y#&PU1D4yb%APKR#!ybk0F09!e=X+IK;SF$AZ0_h@tVjeg^Nc@I66Tvgv$B5MIk
zRkuC)`02AQewyq#^YkaiNA2sFWn;UGPoL;cQSldvn##|#NV9NbbLxZY9y0FJ^
z9X;Sn-1Ex`P29t##ihHla`IS=_sKox{mD@^pIMR1s3=!zL0`BM0}^x0Mry
zMn*a-%=#f}A^#3B(w~F;O5Z(DF7`bwFMir>&ZpR*Ii)L-VsE#f8~r2we6W%OrTXK}
zaG6&GVd2tmXnp!>A2{+cvoLlTD+i^rwDWKhwr*ayNpR}o{yw=-Tad6y(a<&}zoDY8
zz#lp#I%uh-n;8`oO3XwgC--FS>fb`x$Rm!;wv`GI;?7Q;y_bz8b{K)9Vv9>9N1r~G
z)K?3j)#(vYgklH3{ZaH@%V=mAd2RnIDo79R-7;%3NX
zuiYvt48cQTzrC`uvN6+G8+GKiHC^8dX)SCSRWF{+B=9
zKlr17nRDHmp6D*R2k;9(=AiEQz~!(@y69FNy8hTr32z}t5I#yq<73-!ls}Y}T`NAJ
zAR~J_uBY>%xhgsgxsdbgz(6+-zOz1J9zqV)=A8fm$7RLR=0jcZRQvs_ZM`V@E=GEBGlRwn0Io-NEK|_Q`~^25KWJc
z`oiY*>|a5>7iJ6vYI3O(G`zM+_iLs$gZA-MD{B_^s)POtKHW+@Y|^D*iT|l&MJY7I
zM>Xz(voScndHM4v5hkYP(RAS9a@<<$8O>>&dXD%py`Ez@WXmoOnK_S-PxmRcl%lDH
zWq?j3pg@g)EgS?F9VB=p3`#|z1fukLATXbj!Wcrwh-iywIXmcUB@ws7-^%cS+G
zg{;z7I5ZxA@)Q4^cu)Hr&1`Pv-wiVQ_3<&#BL8}sVuY-yElB4Ng=YYmL%6i9Ea51v
z%jqz2hAW@y^tkx$-BbUYK0B3E@Q|L!2-0Sd48y9}gbSxggdK{anwt4!m6eq8YUq28
zr&>>+y2&&7&ydfD+_I;?fB!w=M=mHXE_Mn2>}36ob=EsOg|JWw=lX}4T`(iaG7VBQpTXg+w7*;+=1Y*C-EuW_*o(K;rgygx|(
z@bZSCq8Fq2yn2QGJd2g3hErLo0w2{^l8;Pgd6kLpUUt)T1RFL+v1{P__r0)f{W^>}
z2qo{4^6}A&&v>TCzgYP_m5#kpf19<|eKo%-2R%iU97XLr>dh!!OiXdwPm0FYciPwG
z@(OhKQbZ|UEk*HIYH(8Q=h7IavIv?_qdepps1)8N?#wnKqdR8(a2cjQC*q35j-X?mjhUr!CY_*KEoNYK;a
zU9+PyocHxX<`7jPQ1HY*tL&31#FvuG1TRu>Ww%#n%*hh~?okz?1@F5&r
z;G$)lO224hR-5=xyrIn^x-CxYU1p}buC8w4-mo3_*<`iN*WTWjxLYp#SB6UArRzXO
z?EWdtdCy$fRSkED-3M}-Vkaw
z#>)ehxtsT!>RbQYpqLev#7?X2Aa4t{T7_E!U2*AJ?I)@?w^?IAL4i)+-ya+Dx^7?V
zfqm&?e4e}4qp4qOU@z=FlM&1Fi`X&5$$47z=e%vK=u$(AB
zztJ@|Az^eO$8$9T`~bHscSd4p#ZxvmHq}}WsqVEUSmVbXw>Q7`UZK2wX7A-y^$CYbE%6ehf
zI%iGISW^^5Wo0H&(Z=DDy@A^`tY=VgtwMR~Pi~$W9}mANuA$`A5BLoTe#=vHm#&2!
zHFggBaC@@_2`MQ;(==+61`&aF0BGojJp?qd633&-{>zIqD{E^<142N)=jZ2Ya{ukA
zH#DL(GbMO~gLcn#7Sut75q<6
z%aD6S|Blo@Vk3j}8A8P4%4xD7vm8pNTfaA4z~CcM^s!>{tgBmh)|s7wL|WkfGviC`
z0{X26JYMIb0|R=ad;<68Q4!btOQ|-JUH_!}c*1s6zu}_gQvF^j6ip2y@@vJPDn)Rq
zg-}*jEWwq2^JcGX=-d4K&%Vv?9p5H}k)b_%j0sjl4h}c~IWr!GNYyJT4~xj_s&fME
z#t-UFu66LMLj01kx3`a}cD1=HWXJnUgtIv}!;S$1yjv0>#zM4nBO{NxW9a7W7M&rr
ztTha5OjL@x=xz9%K_P&Ye?RjcM#ttiFZlqY28;xc{d65ci3KGVFd&G9Ns)Q+o7w69
zpSqbXf^9r0HZqgl<*zhKbZb@{7rM3f>Dp#<1cRCPEYuQLkNY^vhQ`O&@`5^hdhYpM
zp1a-(+AwwbAxFXcqL*m2_*nC%+1Kf!z0d5HVa9FA
zb7F~1QxvITBIu2ogo&=?TcWP`PT#4oz9vJa<^|J)CZ=v
zKPr4LP6ZgbxVd8hS=dv!lzi$+EqKUA#z;*d%5cWZ%e%=HHs(o3MMb4kX|9^;lfr9?
z_4Vu58G$30_0eoAs->Qo?RWDR^MQmss&UiPdVm2&vR^@v!^MBk(XTNI+R4~Yn79xsDg`u&RY9(L
zdDZ#A#=avo_d1>ZFXEz)G9fb2#?miduq~8n=+U#GOG$R!Sjx}Ohn@}6&z+T}>v|Be
z^*^TgHEd2SOl-HDarGy1uQngeE{&I)Vlz+*y`)Es^^rI3n*9Ejym{k+mQqAy
z-uQ1Gd@w?e%VA8H*YGLZ;8jGMxhpLOw?;h0ql_tSC#$}C03L+`oO%3Ll;~Fo4s0A8
zP9ogEq#!H^5>8sbC8*J!o{3mFl?#?8H&VpIBM
z{lK-a2$C)=$Hr7GVN8##udjrCZKL**vwYR|WOXDkgD6mSBCF^u7+Cr(w~(Ei+imx1
z(*`~Ww}R5d!VP~dCH#*0NbZ(-X!6T-^0^MY4r_*c0GwAhcNNh|v$dh@x*2-{TCrYe
z$$6-(y#BA{4lLQCm>jlpyrIwi#R!h#CZEQOjDEAni`(W3_(oBmd>tTkK7-up8*z7b
z*6+&4BH1bhpGG8pdV1QP*sTO_|KhN=bp^;Ylf}+(v5Wm-Z@68~n#mJ({9uk;fX0O*xBD@))5E
z7ph;Il(aZ&c(3KS+)I3e+du`%BpGY~G2cvM#Jlg-PNNO6GbIxRK=&RL=Lqbas<8u1
zT`22>Bndf5DOie6CZBPGa|pBXQ7)2zV@P#qnP3t016moNp|gi=%kjeMBTRXii&!l(ZXo~lInB@y%1m;6cB0(c20<=rLU^jmJ
z^yvrx^$s3DO{P%mNF7n4SX*0FDGH
zu5JRL#Kn>D>UQ?eZdF>K_vjJG$lXgj`@O=51a|_tFReFhR%)0Z4#hFFp{?e<
zU?Xb|j(-P22M#HUykP~AZ5?4Jg$bjHB@DY9oG2C{skP8qb67-O+v$OG&?
zk^Y_-I`t<{o+Nf->NI=aM+VpcL<71J%y9$A@#oqxPGNx*406rizkiE9el&r1O-4pm
zMQjKGGW@}KjomE}5!Em{5Vb=|>FBpu_VTMH+C!BFg}MdLUVFlQDK|vTnM8HOq;0A}q#B7h+O_*L>9<{E5
z4C7I@D%z;WRY+s%MmIc5_A<1yf=2oI_gvhf9`-9WEj7cZD!^m^g@X~
znpYX-M|VI716d9+e=J}M<
zkLU?;!^&;*gkF8r~r^cl(%uKjv?Y69pKV_7{PzfaOY^?ir}
z3Aa$gtA&EpsslEO^ouPU%aq>_s`i~Rr2OQSzJed}1WM45vGvCx@xCE6;bwkipC*C%
z2bO1FtPf>q@UzqESefCrm<@3{G6I#t0#a%aa;JG$%+HRh-0H$ZYeVZ5RPp&A7?L<2iA&R7euTT&Ox&O;@
z@kK#Mkl!^6tWT{gH?sPBfNTS_RmI^>ehARDsJ5txvJ$-W|d
zF9iW<@J3;MLx}=_xdOQG8{4*QtRDf-%%gU;G3)Re)VXLn23wI_4MiPJ_wWOOyzFrT
zw}Ud^655y^k2VS!INI{Y#@TyRfJ#Ae{B<{{th6*5h+01Pe?}n$G?V$|8%wr>S9&+m
z5NOK5Av@-R-$4&D888vTPJg)-*-v^TFBhmyJ9sEXoY$UBsU5NuSc$MRa3w~9$Ze4E
z<0RCl#i4%XMuA#-#MI`KE|fF&8q~V7*}?x^YE~rd8G9)U+D$H3<)MR>YzX7vbeN#6HaZa4+^wU2trL|m!OZ-E0fo54fCIX34oo8
z>o*MfV&LU1?>GJ1uqD=cdU44WIOf6aNm1ayYCa_X{pis*xhJUp|`pA<%H4zs*-csZ10No`J0a+es$%
zpTU3Ph#IQwF#oH1JvG#}QX{%tppHSotKVu^ZxCY{D%kc?uT6iaMM4ld0jR*YbuRoT
z8yg#p_53imd^L5S33Lx+C;xDnsgDM&VsbWP_9icW#~;a3E8Uokn3qB-Jd|)-=KjZd
zaNdCEXun-L6i0K5>;iDm&e=utthCH9f0IF31ez9tbu*>KA0IIBt1^+{?tt1N+zvzp
z*e6C8o>;Lu!UT*gxp)`;^Erc~U-uwv*aJ%qjW>eVZK
z-?wpES~@ztFnsMhVE!ZlJKb9M2oPt|6$3`w%UFqI=8w)M8r{mAIA`nybLFELfjXaa
zZ?;7rg}t0MR$%omA)y;w9-bbq7Az#;z#S9nl=L=mxa8kIbU6blcL!b^ma5U^JeocT
z*^!75XnU{wrrec+($`Qr-(ILdl)xfaHqf=S!{ab3M(edv9%ns+#3q>^w~)eDM(Hb*
zp2vWrCOI>Ke!8`7W$}4$((bhccq_7evcDgSAfr{{)Foa8%|I470q@m>IsRF5XXk^S
zag?`laZ8ULZ*#gfej(?0qG4b#09C+LGoeL$shQD{iG!sj#HFL89tCV^Gqc-(XqCbcU!scPeE#LfB$;IOOGj(w+RU$
z1y*3c86Fy%-+#TWwy_6Si<<+)O>d8Y*hypk)xAXn91hvCp=nOc`OEK=x&Bm0e_}1a
zP61rW=QmHqjZr=yU4Xdy`|@mi`vA>D23w3)y=*`_8--eS0jv|JojwbW^^d5O>6l2~
zKl=6H!-q&Bd?hIp#g~83^aGEDqQ30|nJ&$L=ZDMd=QlBt+D-a`@qj5Wz6qQT16SzF
ztE=}Q5(^j~(a>mPe*xgx4toA|9hnBE3<>tz$fJgh{c<