---
title: "process_sleap_data"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{process_sleap_data}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r setup}
library(fishmechr)
library(ggplot2)
library(tidyr)
library(dplyr)
```

```{r}
library(cli)
options(cli.progress_show_after = 0)
options(cli.progress_clear = FALSE)
```

These are your SLEAP output data files. This is example code that you could use. For the vignette, these data files are included in the package.
```{r, eval=FALSE}
library(here)

sleapfiles <- c(
  here("data-raw","2024-11-14_labels2.000_fish_01-030RPM-ortho-2024-11-14T044441.analysis.csv"),
  here("data-raw","2024-11-14_labels2.001_fish_01-040RPM-ortho-2024-11-14T044525.analysis.csv"),
  here("data-raw", "2024-11-14_labels2.002_fish_01-050RPM-ortho-2024-11-14T045204.analysis.csv"),
  here("data-raw","2024-11-14_labels2.003_fish_01-060RPM-ortho-2024-11-14T044754.analysis.csv"),
  here("data-raw","2024-11-14_labels2.004_fish_01-070RPM-ortho-2024-11-14T044858.analysis.csv")
)

parse_file_name <- function(fn)
{
  d <- stringr::str_match(fn, 
                 "_(?<id>\\w+_\\d+)-(?<speed>\\d+)RPM.+(?<datetime>\\d{4}-\\d{2}-\\d{2}T\\d{6})")
  
  tibble::as_tibble_row(d[1, 2:4], .name_repair = "minimal")
}

zfishdata <- purrr::map(sleapfiles, \(fn) readr::read_csv(fn, id = "fn", 
                                         show_col_types = FALSE)) |> 
  bind_rows() |> 
  mutate(dd = purrr::map(fn, parse_file_name)) |> 
  unnest(dd) |> 
  mutate(fn = basename(fn),
         speed = as.numeric(speed)) 

zfish_goodframes <- readr::read_csv(here("data-raw", "zfish_goodframes.csv"))
```

```{r}
head(zfishdata)
head(zfish_goodframes)
```


Frame rate from your cameras
```{r}
fps <- 50
```

These are all our points, ordered from head to tail
```{r}
pointnames <- c("snout", "eye_ctr", "left_eye", "right_eye", "hyoid",   
                "pec_fin_ctr", "left_pec_fin_base", "left_pec_fin_tip", 
                "right_pec_fin_base", "right_pec_fin_tip", "pelvic_fin_base",
                "anus", "peduncle", "ventral_caudal", "dorsal_caudal")
```

```{r}
zfdata_ctr <-
  zfishdata |>  
  mutate(eye_ctr.x = (left_eye.x + right_eye.x)/2,
         eye_ctr.y = (left_eye.y + right_eye.y)/2,
         eye_ctr.score = NA,
         pec_fin_ctr.x = (left_pec_fin_base.x + right_pec_fin_base.x)/2,
         pec_fin_ctr.y = (left_pec_fin_base.y + right_pec_fin_base.y)/2,
         pec_fin_ctr.score = NA
         )
```


```{r}
zfdata_ctr <- zfdata_ctr |> 
  relocate(id, speed, datetime) |> 
  pivot_kinematics_longer(pointnames = pointnames)
```

```{r}
head(zfdata_ctr)
```

Now pull out just the good frames based on our other data table
```{r}
zfdata_good <- list()
for (i in seq(1, nrow(zfish_goodframes))) {
  good1 <- zfdata_ctr |> 
    filter(fn == zfish_goodframes$File[[i]],
           between(frame_idx, zfish_goodframes$Start[[i]], zfish_goodframes$End[[i]])) |> 
    mutate(block = zfish_goodframes$Block[[i]])
  zfdata_good[[i]] <- good1
}
zfdata_good <- bind_rows(zfdata_good)
```

This shows the x and y positions of several of the points in multiple frames.
```{r}
zfdata_good |> 
  arrange(speed, frame_idx, point) |> 
  filter(speed == 30,
         between(frame_idx, 120, 140)) |> 
  filter(point %in% c("snout", "hyoid", "eye_ctr", "pec_fin_ctr", 
                      "pelvic_fin_base", "anus", "peduncle", "ventral_caudal")) |> 
  ggplot(aes(x = x, y = y, color = point)) +
  geom_point() +
  geom_path(aes(group = frame_idx)) +
  coord_fixed()
```

This shows just the y coordinate against frame number.
```{r}
zfdata_good |> 
  arrange(speed, frame_idx, point) |> 
  filter(speed == 40,
         between(frame_idx, 120, 250)) |> 
  filter(point %in% c("snout", "peduncle", "ventral_caudal")) |> 
  ggplot(aes(x = frame_idx, y = y, color = point)) +
  geom_point() +
  geom_path()
```

This smooths the data and fills in gaps. Focus just on the points along the ventral midline for right now.
```{r}
zfdata_sm <-
  zfdata_good |> 
  arrange(speed, frame_idx, point) |> 
  filter(point %in% c("snout", "hyoid", "eye_ctr", "pec_fin_ctr", 
                      "pelvic_fin_base", "anus", "peduncle", "ventral_caudal")) |> 
  group_by(fn, frame_idx) |>
  # calculate the arc length
  mutate(arclen0 = arclength(x, y, na.skip = TRUE)) |> 
  # smooth and fill gaps
  interpolate_points_df(arclen0, x, y, spar = 0.2,
                        tailmethod = 'extrapolate',
                        fill_gaps = 1,
                        .frame = frame_idx,
                        .out = c(arclen='arclen', xs='xs', ys='ys')) |> 
  ungroup()
  
```

```{r}
zfdata_sm |> 
  arrange(speed, frame_idx, point) |> 
  filter(speed == 40,
         between(frame_idx, 120, 250)) |> 
  filter(point %in% c("snout", "peduncle", "ventral_caudal")) |> 
  ggplot(aes(x = frame_idx, y = ys, color = point)) +
  geom_point() +
  geom_path()

```

Now we need to find the center of mass. For that we need the width and height of the body at each point along the body. These values are all in fractions of the total length.
```{r}
zebrafish_shape |> 
  ggplot(aes(x = s)) +
  geom_path(aes(y = width)) +
  geom_path(aes(y = height), color = "blue")
```

Now interpolate and scale these width and height values for the actual length of our fish.
```{r}
zfdata_sm <- zfdata_sm |> 
  group_by(id, speed, datetime, frame_idx) |> 
  mutate(width = interpolate_width(zebrafish_shape$s,
                                   zebrafish_shape$width,
                                   arclen),
         height = interpolate_width(zebrafish_shape$s,
                                   zebrafish_shape$height,
                                   arclen)
         )
```

```{r}
zfdata_sm |> 
  filter(speed ==30,
         frame_idx == 150)
```

Now we need to split out different swimming sequences
```{r}
zfdata_split <-
  zfdata_sm |> 
  group_by(id, speed, datetime, block) |> 
  group_split()
```

```{r}
zfdata_ctr <- list()
for (i in seq(1, length(zfdata_split))) {
  zfdata_ctr[[i]] <-
    zfdata_split[[i]] |> 
    get_midline_center_df(arclen, xs, ys,
                          width = width, height = height,
                          .frame = frame_idx) |> 
    # center everything on the center of mass
    mutate(xctr = xs - xcom,
           yctr = ys - ycom,
           t = frame_idx / fps) |> 
    # find the main axis of the body
    get_primary_swimming_axis_df(t, xctr, yctr, 
                                 .frame = frame_idx)
}
```

```{r}
zfdata_ctr[[2]] |> 
  ggplot(aes(x = exc_x, y = exc, color = point)) +
  geom_path(aes(group = frame_idx)) +
  geom_point()
```

```{r}
zfdata_ctr[[2]] |> 
  filter(point %in% c("peduncle", "ventral_caudal")) |> 
  # filter(between(t, 1, 1.2)) |> 
  ggplot(aes(x = t, y = exc, color = point)) +
  geom_path()
```

```{r}
zfdata_phase <- list()

for (i in seq(1, length(zfdata_ctr))) {
  zfdata_phase[[i]] <- 
  zfdata_ctr[[i]] |> 
  arrange(id, speed, datetime, frame_idx, desc(point)) |> 
  group_by(id, speed, datetime, point) |> 
  mutate(ph_p = peak_phase(exc))
}
```

```{r}
zfdata_phase[[2]] |> 
  ungroup() |> 
  filter(point %in% c("snout", "peduncle", "ventral_caudal")) |> 
  mutate(point = factor(point)) |> 
  ggplot(aes(x = t, y = ph_p, color = point)) +
  geom_path() + 
  facet_wrap(~point)
```

```{r}
zfdata_phase[[2]] |> 
  ungroup() |> 
  group_by(point) |> 
  mutate(freq_p = get_frequency(t, ph_p, method='deriv')) |> 
  filter(point %in% c("snout", "peduncle", "ventral_caudal")) |> 
  ggplot(aes(x = t, y = freq_p, color = point)) +
  scale_shape_manual(values = c(1, 17, 22)) +
  geom_point() +
  facet_wrap(~point)
```

```{r}
zfdata_cyc <- list()

for (i in seq(1, length(zfdata_phase))) {
  zfdata_cyc[[i]] <- 
    zfdata_phase[[i]] |> 
    group_by(point) |> 
    mutate(freq = get_frequency(t, ph_p, method='deriv')) |>
    ungroup() |> 
    get_body_cycle_numbers_df(ph_p, pointval = "peduncle",
                              .frame = frame_idx) |> 
    arrange(id, speed, datetime, block, t, point)
}
```

```{r}
zfdata_cyc[[1]] |> 
  ungroup() |> 
  group_by(speed, point, cycle) |> 
  summarize(amp = (max(exc) - min(exc)) / 2,
            arclen = mean(arclen)) |> 
  ggplot(aes(x = arclen, y = amp, color = speed)) +
  geom_path(aes(group = cycle))

```

```{r}
zfdata_cyc <- bind_rows(zfdata_cyc)
```

```{r}
zfsummary <-
  zfdata_cyc |> 
  group_by(id, datetime, speed, frame_idx) |> 
  mutate(meanfreq = mean(freq)) |> 
  filter(point == "ventral_caudal") |> 
  group_by(speed, cycle) |> 
  summarize(amp = (max(exc) - min(exc)) / 2,
            arclen = mean(arclen),
            meanfreq = mean(meanfreq)) |> 
  ungroup()

zfsummary |> 
  ggplot(aes(x = speed, y = meanfreq)) +
  geom_point(aes(group = factor(speed)))
```