---
title: "Diaptera wings classification using Topological Data Analysis"
author:
  - name: Guilherme Vituri F. Pinto
    orcid: 0000-0002-7813-8777
    corresponding: true
    email: vituri.vituri@gmail.com
    roles:
      - Investigation
      - Project administration
      - Software
      - Visualization
    affiliations:
      - Unesp
  - name: Sergio Ura
  - name: Northon
keywords:
  - Topological Data Analysis
  - Persistent homology
abstract: |
  We studied etc etc etc  etc etc etc  etc etc etc  etc etc etc  etc etc etc  etc etc etc  etc etc etc  etc etc etc  etc etc etc  etc etc etc  etc etc etc  etc etc etc  etc etc etc  etc etc etc 
plain-language-summary: |
  Earthquake data for the island of La Palma from the September 2021 eruption is found ...
key-points:
  - A web scraping script was developed to pull data from the Instituto Geogràphico Nacional into a machine-readable form for analysis
  - Earthquake events on La Palma are consistent with the presence of both mantle and crustal reservoirs.
date: last-modified
bibliography: references.bib
citation:
  container-title: Earth and Space Science
number-sections: true
---

```{julia}
using TDAfly, TDAfly.Preprocessing, TDAfly.TDA, TDAfly.Analysis
using Images: mosaicview
using Plots: plot, display, heatmap, scatter
using PersistenceDiagrams
```

## Introduction

Falar sobre o dataset, TDA, etc.

## Methods

All images are in the `images/processed` directory. For each image, we load it, apply a gaussian blur, crop and make it have 150 pixels of height. The blurring step is necessary to "glue" small holes in the figure and keep it connected.

```{julia}
paths = readdir("images/processed", join = true)
species = basename.(paths) .|> (x -> replace(x, ".png" => ""))
individuals = map(species) do specie
  s = split(specie, " ")
  s[1][1] * "-" * s[2]
end
wings = load_wing.(paths, blur = 1.3)
Xs = map(image_to_r2, wings);
```

```{julia}
mosaicview(wings, ncol = 4, fillvalue=1)
```

### Vietoris-Rips filtration

We select 500 points from each image using a farthest point sample method

```{julia}
samples = map(Xs) do X
  ids = farthest_points_sample(X, 500)
  X[ids]
end;
```

We then calculate its persistence diagrams using the Vietoris-Rips filtration etc.

```{julia}
# get only the 1-dimensional PD
rips = rips_pd.(samples, cutoff = 5, threshold = 200) .|> last;
```

We create the 1-dimensional persistence image for each persistence diagram using 10x10 matrices

```{julia}
PI = PersistenceImage(rips, size = (10, 10))

images_rips = PI.(rips);
```

### Examples

Below are some examples of 1-dimensional barcodes, its persistence image and the original wing that generated it. Note: we are plotting the barcode using the birth and persistence.

```{julia}
# plot some images to see the barcodes
map([1, 4, 8, 10, 15]) do i
  p = plot_wing_with_pd(rips[i], images_rips[i], samples[i], species[i])
  display(p)
end;
```


We now calculate the Euclidean distance between each persistence image (seen as a vector of $\mathbb{R}^{10x10}$) and plot its heatmap

```{julia}
D_rips = pairwise_distance(images_rips);
```

```{julia}
plot_heatmap(D_rips, individuals, "Distance matrix for Vietoris-Rips barcodes")
```

### Persistence Homology Transform

Now we will create several filtrations based on points and lines, etc.

We start with the point (0, 0). Its filtration is the following

```{julia}
A = wings[1] |> image_to_array;
f = dist_to_point(0, 0)
Af = modify_array(A, f)
heatmap(Af)
```

with corresponding sublevel barcode as

```{julia}
point_pds = cubical_pd(Af, cutoff = 0.05)
plot_pd(point_pds)
```

or, with persistence in the y-axis:

```{julia}
plot_pd(point_pds, persistence = true)
```

Let's see step-by-step of this filtration:

```{julia}
for tr ∈ reverse([0:0.1:1;])
  X = findall_ids(>(tr), Af)
  title = "threshold: $tr"
  p = scatter(first.(X), last.(X), title = title)
  display(p)
end
```

Due to noise, some connected components are born in 0.2 and die only at 0. But the loops seems alright.



## Draft.....

```{julia}
#| include: false 
ps = map(reverse([0:0.1:0.9;])) do tr
  X = findall_ids(>(tr), Af)
  title = "threshold: $tr"
  p = scatter(first.(X), last.(X), title = title)
  # display(p)
end

plot(ps...)
```