04bis.matches_stats_apeep_manual.Rmd

---
title: "Statistics for segmentation benchmark (Apeep VS manual on CC4)"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.width=10, fig.height=6)
library(tidyverse)
library(scales)
```

The purpose of this script is to compare the performances of apeep segmentation against a ground truth human segmentation, on 258 large `apeep` images (10240 px * 2048 px).  
Apeep pipeline is an adaptive gray level segmentation, down to 50 px.  


## Prepare data
### Read data


```{r read_data}
output_dir <- "data_cc4/matches"

man_parts <- read_csv(file.path(output_dir, "man_particles_props.csv"), col_types = cols()) %>% select(-c("object_label", "object_bbox-0", "object_bbox-1", "object_bbox-2", "object_bbox-3"))
reg_parts <- read_csv(file.path(output_dir, "reg_particles_props.csv"), col_types = cols()) %>% select(-c("object_label", "object_bbox-0", "object_bbox-1", "object_bbox-2", "object_bbox-3"))

matches_reg <- read_csv(file.path(output_dir, "matches_reg.csv"), col_types = cols())
```


### Select relevant objects

Make a list of taxa in manually segmented particles.

```{r taxa}
taxa <- man_parts %>% pull(taxon) %>% unique() %>% sort()
taxa
```

The manual segmentation originally generated `r nrow(man_parts)` particles.
We will ignore objects in the `detritus` and `othertocheck` categories as well as objects smaller than 50 px.

```{r filter_objects}
ignored <- man_parts %>% filter(taxon %in% c("detritus", "othertocheck")) %>% pull(object_id)
small <- man_parts %>% filter(area < 50) %>% pull(object_id)
man_parts <- man_parts %>% filter(!(taxon %in% c("detritus", "othertocheck"))) %>% filter(area >= 50)

matches_reg <- matches_reg %>% filter((!(man_ids %in% ignored)) & (!(man_ids %in% small))) 
```

**After removing non living and small particles, `r nrow(man_parts)` manual particles are left.
`r nrow(reg_parts)` particles were generated by apeep segmentation.**

Let’s inspect taxonomic composition of benchmark dataset.

```{r testset_comp}
man_parts %>% 
  count(taxon) %>% 
  arrange(-n) %>% 
  mutate(taxon = factor(taxon, taxon)) %>% 
  ggplot() +
  geom_col(aes(x = taxon, y = n, fill = n > 10)) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_y_continuous(trans = "log1p", breaks = c(0, 10, 50, 100, 200, 400, 1000, 2000)) +
  labs(x = "Taxon", y = "Object number", title = "Test set composition") + 
  theme(text = element_text(size = 16)) 
  
```


## Compute global statistics
Compute overall precision and recall:  

* precision: among `apeep` particles, how many where matched with manual particles?
* recall: among manual particles, how many where matched with `apeep` particles?

```{r global}
# Precision
# apeep particles matched with manual particles / apeep particles
precision <- length(unique(matches_reg$reg_ids)) / length(reg_parts$object_id)

# Recall
# manual particles matched with apeep particles / manual particles
recall <- length(unique(matches_reg$man_ids)) / length(man_parts$object_id)

global_stats <- tibble(
  metric = c("precision", "recall"),
  value = c(precision, recall)
)
global_stats
```

**Apeep segmentation has a very good recall (`r percent(recall, accuracy = 0.1)`) but a very poor precision (`r percent(precision, accuracy = 0.1)`).**


## Compute statistics per taxon
We want to compute the recall of organism per taxonomic group, but we have to deal with multiple matches.

Case 1: one `apeep` particle matched with multiple manual particles, likely with two different taxo.
Two solutions:

- take the rarest taxo
- ignore particle as the CNN won’t be able to predict it (selected solution)

Case 2: one manual particle matched with multiple `apeep` particles, only one taxo but `apeep` segmentation overestimates the number of organisms in this taxo. 

Solution: keep only one match.

```{r taxo_stats}
## apeep particles
# Count matches 
counts_taxo <- matches_reg %>% 
  # drop all cases of duplicated matches of apeep particles in match table (solve case 1 of multiple matches)
  add_count(reg_ids) %>% filter(n==1) %>% select(-n) %>% 
  # drop duplicates of matched manual particles and keep only one (solve case 2 of multiple matches)
  distinct(man_ids, .keep_all = TRUE) %>% 
  # join matched apeep particles with manual particles taxo
  left_join(man_parts, by = c("man_ids" = "object_id")) %>% 
  # count matches per taxon
  count(taxon, name = "n_reg") 

# Compute recall
recall_taxo <- man_parts %>% 
  # count true particles per taxon
  count(taxon, name = "n_truth") %>% 
  # join with matched particles
  left_join(counts_taxo, by = "taxon") %>% 
  # compute ratio of apeep matched particles over true particles (recall)
  mutate(recall= n_reg / n_truth) %>% 
  arrange(n_truth) %>% 
  mutate(taxon = factor(taxon, taxon))

# Plot it
recall_taxo %>% 
  ggplot() +
  geom_col(aes(x = taxon, y = n_truth, fill = recall), position = "dodge") +
  scale_fill_viridis_c(direction = -1) + 
  scale_y_continuous(trans = "log1p") +
  coord_flip() +
  labs(title = "Recall scores and organisms number per taxa for each segmentation", y = "Number of true objects") +
  theme(text = element_text(size = 16)) 
```

Other plot 

```{r taxo_plot, echo=FALSE}
recall_taxo %>% 
  ggplot() +
  geom_col(aes(x = taxon, y = recall, fill = n_truth)) + 
  scale_fill_viridis_c(trans = "log1p", direction = -1) + 
  coord_flip() +
  labs(title = "Recall scores and organisms number per taxa", fill = "Number of \ntrue objects") +
  theme(text = element_text(size = 16)) 

```

## Compute statistics per size class

We will define size classes for particles:

* [50 px,100 px)  
* [100 px, 150 px) 
* [150 px, 200 px) 
* [200 px, 250 px) 
* [250 px, 300 px) 
* [300 px, 350 px) 
* [350 px, 400 px) 
* [400 px, 450 px) 
* [450 px, 500 px) 
* \> 500 px

And compute recall for each size class.


```{r size_stats}
# Define size classes
# - 50 px from 50 px to 500 px
# - larger than 500 px
man_parts <- man_parts %>% 
  mutate(class_size = cut(area, breaks = c(0, seq(from = 50, to = 500, by = 50), 1000000), right = FALSE))

## apeep particles
# Count matches 
counts_size <- matches_reg %>% 
  # drop all cases of duplicated matches of apeep particles in match table (solve case 1 of multiple matches)
  add_count(reg_ids) %>% filter(n==1) %>% select(-n) %>% 
  # drop duplicates of matched manual particles and keep only one (solve case 2 of multiple matches)
  distinct(man_ids, .keep_all = TRUE) %>% 
  # join matched apeep particles with manual particles class size
  left_join(man_parts, by = c("man_ids" = "object_id")) %>% 
   # count matches per class size
  count(class_size, name = "n_reg")

# Compute recall
recall_size <- man_parts %>% 
  # count true particles per class size
  count(class_size, name = "n_truth") %>% 
  # join with matched particles
  left_join(counts_size, by = "class_size") %>% 
  # compute ratio of apeep matched particles over true particles (recall)
  mutate(recall = n_reg / n_truth) 

recall_size %>% 
  ggplot() +
  geom_col(aes(x = class_size, y = n_truth, fill = recall), position = "dodge") +
  scale_fill_viridis_c(direction = -1) + 
  scale_y_continuous(trans = "log1p") +
  coord_flip() +
  labs(title = "Recall scores and organisms number per size class", fill = "Recall", 
       x = "Size class (px)", y = "Number of \ntrue objects") +
  theme(text = element_text(size = 16)) 
```

Other plot.

```{r size_plot, echo=FALSE}
recall_size %>% 
  ggplot() +
  geom_col(aes(x = class_size, y = recall, fill = n_truth)) + 
  scale_fill_viridis_c(trans = "log1p", direction = -1) + 
  coord_flip() +
  labs(title = "Recall scores and organisms number per size class", fill = "Number of \ntrue objects", x = "Size class (px)") +
  theme(text = element_text(size = 16)) 

```



Now compute precision on each size class.

```{r size_prec}
# Define size classes
# - 50 px from 50 px to 500 px
# - larger than 500 px
reg_parts <- reg_parts %>% mutate(class_size = cut(object_area, breaks = c(0, seq(from = 50, to = 500, by = 50), 1000000), right = FALSE))

## apeep particles
# Count matches 
counts_size <- matches_reg %>% 
  # drop all cases of duplicated matches of apeep particles in match table (solve case 1 of multiple matches)
  add_count(reg_ids) %>% filter(n==1) %>% select(-n) %>% 
  # drop duplicates of matched manual particles and keep only one (solve case 2 of multiple matches)
  distinct(man_ids, .keep_all = TRUE) %>% 
  # join matched apeep particles with apeep particles class size
  left_join(reg_parts, by = c("reg_ids" = "object_id")) %>% 
   # count matches per class size
  count(class_size, name = "n_match")

# Compute precision
precision_size <- reg_parts %>% 
  # count true particles per class size
  count(class_size, name = "n_reg") %>% 
  # join with matched particles
  left_join(counts_size, by = "class_size") %>% 
  # compute ratio of apeep matched particles over true particles (recall)
  mutate(precision = n_match / n_reg) %>% 
  # add count of true objects
  left_join(man_parts %>% count(class_size, name = "n_truth"),  by = "class_size") 

precision_size %>% 
  ggplot() +
  geom_col(aes(x = class_size, y = n_truth, fill = precision), position = "dodge") +
  scale_fill_viridis_c(direction = -1) + 
  scale_y_continuous(trans = "log1p") +
  coord_flip() +
  labs(title = "Precision scores and organisms number per size class for each segmentation", fill = "Precision", 
       x = "Size class (px)", y = "Number of \ntrue objects") +
  theme(text = element_text(size = 16)) 
```


Other plot.

```{r size_prec_plot, echo=FALSE}
precision_size %>% 
  ggplot() +
  geom_col(aes(x = class_size, y = precision, fill = n_truth)) + 
  scale_fill_viridis_c(trans = "log1p", direction = -1) + 
  coord_flip() +
  labs(title = "Precision scores and organisms number per size class", fill = "Number of \ntrue objects", x = "Size class (px)") +
  theme(text = element_text(size = 16)) 

```