We calculated the reflectance of the beetles from the RGB values measured from calibrated photographs to find the degree of circular polarisation.

First we linearised the RBGs using the parameters derived from the MATLAB code adapted from Stevens et al 2007 (Biological Journal of the Linnean society, 90(2), 211-237). Next we equalised the RGB values with the grey standard with known reflectance. Finally, we calculated the difference between the two Polarisations standardized by their sum.

In the second part of this script, we correlate the polarisation with reflectivity obtained in previous scripts

Setting up

General

Libraries sourced from an additional script

source(here::here("Scripts/MacroEcol_1_Libraries.R"))

We included an overall option to keep as many decimals as possible in our data, which is important for our calculations

options(scipen = 100)

We imported the data set containing ecological and phylogeny data obtained in previous steps.

Cons1 <- read.csv(here::here("Data/FromCode/ConsReflEcolInd.csv"))[-1]

Parameters

To calibrate, we photographed a set of 2%, 20%, 40%, 60%, 80% and 99% and standards. Then, we fitted a polynomial regression between the RGB values the standards display in our set up and their expected reflectance on each channel (obtained from the manufacturer), using a custom made function in MatLab.

Since the pictures were not taken on the same day, we calibrated the beetle photographs to the photographs of the standards taken on the same day. That means we had one set of parameters for each camera configuration (in the code this is called camera_cat)

The parameters obtained from these curves were stored in .txt files, so we wrote a function to import them directly.

read_parameters_from_file <- function(camera_cat, type) {
  param <- read.table(paste0(
    here::here("Data/ParametersMatLab/StandardParams"),
    camera_cat, type, ".txt"
  ), header = TRUE) %>%
    dplyr::select(1:5, 7) %>%
    dplyr::mutate(camera_cat = rep(camera_cat)) %>%
    dplyr::mutate(type = rep(type))

  return(param)
}

parameter0 <-
  bind_rows(
    read_parameters_from_file("J21", "LCP"),
    read_parameters_from_file("J21", "RCP"),
    read_parameters_from_file("J21", "VIS"),
    read_parameters_from_file("J25", "LCP"),
    read_parameters_from_file("J25", "RCP"),
    read_parameters_from_file("J25", "VIS"),
    read_parameters_from_file("J26", "LCP"),
    read_parameters_from_file("J26", "RCP"),
    read_parameters_from_file("J26", "VIS"),
    read_parameters_from_file("F01", "LCP"),
    read_parameters_from_file("F01", "RCP"),
    read_parameters_from_file("F01", "VIS"),
    read_parameters_from_file("F04", "LCP"),
    read_parameters_from_file("F04", "RCP"),
    read_parameters_from_file("F04", "VIS")
  )

The regression uses 4 parameters a, b, c, and d. We obtained one value for each parameter, on each channel (R,G, B) for each type of filter (VIS = Visible, LCP =left handed circularly polarized and RCP= right circularly polarized) on each camera configuration (camera_cat)

parameter <-
  parameter0 %>% dplyr::filter(Channel != "Gy") %>% # keep only the channels R, G, B.
  dplyr::select(-5) %>% # this column is the % of fitting (ideally) above 0.99)
  gather(
    key = parameter,
    value = value,
    -Channel, -type, -camera_cat
  ) # re-arrange: All parameters in 1 column

head(parameter)

Channel	camera_cat	type	parameter	value
R	J21	LCP	a	2.38e+03
G	J21	LCP	a	0.000544
B	J21	LCP	a	0.0275
R	J21	RCP	a	0.0105
G	J21	RCP	a	6.68e-06
B	J21	RCP	a	0.000203

Grey Standard

This is the known reflectance (from manufacturer) of the 40% grey placed next to each beetle in the photographs

grey_standard_refl <- read.csv(here::here("Data/12_knownGreyStdReflectance.csv"))

Raw RGB values

These are the RGB valuyes extracted from the calibrated photographs. We considered one region of interest (ROI) in one elytron and one in the pronotum. We also sampled the 40% grey as a reference.

Elytra

We separated the RGB values of the elytra and the grey standard but kept a unique identifier called PhotoID. This helps pairing the sample with the grey that was in that picture further in the analysis (required to run the linearisation function).

RGBElytraNGrey <-
  RGB_raw_transform %>%
  dplyr::filter(tr != "Pronotum") %>% # remove data for the pronotum
  unite("PhotoID", # Create a unique code
    c(1, 2, 4), # which contains the individual_Filter_Channel
    remove = FALSE
  ) # still keep all the original columns

# Each beetle sample has a grey standard in the same photo.
levels(as.factor(RGBElytraNGrey$tr))

## [1] "Elytra" "Grey"

# We will separate the elytra and grey
# but they will still share the unique PhotoID to pair them later

RGBElytra <-
  RGBElytraNGrey %>%
  dplyr::filter(tr == "Elytra")

RGBElGrey <-
  RGBElytraNGrey %>%
  dplyr::filter(tr == "Grey")

Pronotum

Pronotum data was not available for the following beetles:

RGBPronot1 <-
  RGB_raw_transform %>%
  dplyr::filter(tr != "Elytra") %>% # remove data for the elytra
  unite("PhotoID", # Create a unique code
    c(1, 2, 4), # which contains the individual_Filter_Channel
    remove = FALSE
  ) # still keep all the original columns

inPGind <- levels(as.factor(RGBPronot1[RGBPronot1$tr == "Grey", ]$ind))
inPRind <- levels(as.factor(RGBPronot1[RGBPronot1$tr == "Pronotum", ]$ind))

setdiff(inPGind, inPRind)

## [1] "ecry03" "mcla02" "roci02"

setdiff(inPRind, inPGind)

## character(0)

# Note pronotum missing for these: "ecry03" "mcla02" "roci02"

They are considered for the overall polarisation analysis but not in the correlations elytra-pronotum.

We separated the RGB values of the pronotum and the grey standard but kept a unique identifier called PhotoID. This helps pairing the sample with the grey that was in that picture further in the analysis (required to run the linearisation function).

RGBPronot1 <-
  RGBPronot1 %>%
  dplyr::filter(ind != "ecry03" &
    ind != "mcla02" &
    ind != "roci02") # Remove these species because pronotum unavailable

RGBPronot <-
  RGBPronot1 %>%
  dplyr::filter(tr == "Pronotum")

RGBPrGrey <-
  RGBPronot1 %>%
  dplyr::filter(tr == "Grey")

Calibrations

Here we linearise and equalise RGBs for each camera setting and each channel separately.

First we defined the functions as follows:

Linearization

\[\ Y = a*exp(b*x) + c*exp(d*x)\] Where a, b c and d are the parameters obtained from MatLab. x is the raw RGB value for the ROI and Y is the linearised value

# The arguments of this functions are:
# a DataFrame: containing the RGB values that we want to linearise
# a ParameterData: a data frame containing the parameters a,b,c,d from MatLab

lin_equ_RGB <- function(DataFrame, ParameterData) {
  Result <- DataFrame # keep the original data, results will be added to it
  Result$Linearized <- rep("NA", length(DataFrame[1])) # add an empty column

  for (i in 1:length(DataFrame$PhotoID)) {
    para <- ParameterData %>% # the data frame with the parameters

      # For each line in the DataFrame to linearise, we need to find
      # the correspondent set of parameters:
      dplyr::filter(type == DataFrame$type[i], # VIS, LCP, or RCP
        Channel == DataFrame$Channel[i], # R, G or B
        camera_cat == DataFrame$camera_cat[i]) %>% # Camera conditions
      dplyr::arrange(., parameter) # arranged alphabetically

    # Apply the formula
    Result$Linearized[i] <-
      para$value[1] * exp(para$value[2] * Result$value[i]) +
      para$value[3] * exp(para$value[4] * Result$value[i])
  }
  return(Result)
}

Equalisation

\[\ E = l / (rg /(kg*255))\] Where rg is the raw RGB values for the grey in the photo, kg is the known reflectance of the gray standard (from the fabricator), l is the linearised value of the ROI and E is the equalised value of the ROI.

For this part we did not use and equation, rather we obtained the calculations directly form the data frames.

Elytra

Linearise

We applied the function defined in the previous step

LinearizedElytra <- lin_equ_RGB(RGBElytra, parameter) %>%
  dplyr::select(PhotoID, tr, Linearized) # apply function

head(LinearizedElytra)

PhotoID	tr	Linearized
abno01_VIS_B	Elytra	0.0433459224340714
abno01_LCP_B	Elytra	0.06478427380179
abno01_RCP_B	Elytra	0.0491076172571187
abno01_VIS_G	Elytra	0.0739528114682013
abno01_LCP_G	Elytra	0.0978269786698409
abno01_RCP_G	Elytra	0.0954403321081553

LinearizedElytra$Linearized <- round(as.numeric
(LinearizedElytra$Linearized), 4) # round

LinearizedElytraSpr <-
  LinearizedElytra %>%
  spread(tr, Linearized)

Equalise

And then equalised following the equation

RGBElGreys <-
  RGBElGrey %>%
  dplyr::select(PhotoID, value) # select the relevant columns of the raw grey values

EQElytra <-
  merge(RGBElGreys, LinearizedElytraSpr,
    all = FALSE
  ) %>%
  dplyr::rename(RawGrey = value) %>%
  dplyr::mutate(
    Equalized =
      Elytra / (RawGrey / (grey_standard_refl[1, 2] * 255))
  )

Average

We averaged the RGB values of the three channels for each filter.

RGB_ave <- EQElytra %>%
  dplyr::select(PhotoID, Equalized) %>%
  dplyr::mutate(Reflectance = Equalized * 100) %>%
  dplyr::select(-Equalized) %>%
  separate(PhotoID, c("ind", "Type", "Channel"), sep = "_") %>%
  dplyr::group_by(ind, Type) %>%
  dplyr::mutate(RGB_ave = mean(Reflectance)) %>%
  dplyr::select(ind, Type, RGB_ave) %>%
  dplyr::distinct() %>% # leave only the one row for each individual
  dplyr::filter(!grepl("ambl", ind)) %>%
  dplyr::filter(!grepl("sqzb", ind)) # remove, absent in the phylogeny

head(RGB_ave)

## # A tibble: 6 × 3
## # Groups:   ind, Type [6]
##   ind    Type  RGB_ave
##   <chr>  <chr>   <dbl>
## 1 abno01 LCP      7.40
## 2 abno01 RCP      6.61
## 3 abno01 VIS      5.31
## 4 abno02 LCP      5.95
## 5 abno02 RCP      5.15
## 6 abno02 VIS      4.15

Pronotum

Linearise

We applied the function defined in the previous step

LinearizedPron <- lin_equ_RGB(RGBPronot, parameter) %>%
  select(PhotoID, tr, Linearized) # apply function

head(LinearizedPron)

PhotoID	tr	Linearized
rept01_VIS_R	Pronotum	0.0486138501955103
rept01_VIS_G	Pronotum	0.0493577182264623
rept01_VIS_B	Pronotum	0.039301457688771
rept01_LCP_R	Pronotum	0.0519620880912811
rept01_LCP_G	Pronotum	0.0636280773901934
rept01_LCP_B	Pronotum	0.0439852217645462

LinearizedPron$Linearized <- round(as.numeric
(LinearizedPron$Linearized), 4) # round

LinearizedPronSpr <-
  LinearizedPron %>%
  spread(tr, Linearized)

Equalise

And then equalised following the equation

RGBPrGreys <-
  RGBPrGrey %>%
  dplyr::select(PhotoID, value) # select the relevant columns of the raw grey values

EQPronot <-
  merge(RGBPrGreys, LinearizedPronSpr,
    all = FALSE
  ) %>%
  dplyr::rename(RawGrey = value) %>%
  dplyr::mutate(
    Equalized =
      Pronotum / (RawGrey / (grey_standard_refl[1, 2] * 255))
  )

Average

We averaged the RGB values of the three channels for each filter.

RGB_avePr <- EQPronot %>%
  dplyr::select(PhotoID, Equalized) %>%
  dplyr::mutate(Reflectance = Equalized * 100) %>%
  dplyr::select(-Equalized) %>%
  separate(PhotoID, c("ind", "Type", "Channel"), sep = "_") %>%
  dplyr::group_by(ind, Type) %>%
  dplyr::mutate(RGB_aveP = mean(Reflectance)) %>%
  dplyr::select(ind, Type, RGB_aveP) %>%
  dplyr::distinct() %>% # leave only the one row for each individual
  dplyr::filter(!grepl("ambl", ind)) %>%
  dplyr::filter(!grepl("sqzb", ind)) # remove, absent in the phylogeny

head(RGB_avePr)

## # A tibble: 6 × 3
## # Groups:   ind, Type [6]
##   ind    Type  RGB_aveP
##   <chr>  <chr>    <dbl>
## 1 abno01 LCP       6.60
## 2 abno01 RCP       5.68
## 3 abno01 VIS       4.68
## 4 abno02 LCP       5.76
## 5 abno02 RCP       4.81
## 6 abno02 VIS       4.01

Polarisation

We calculated the polarisation as \[\ P \ = \frac{l - r}{l + r}\]

where \(l\) represents the mean reflectance under the filter for left handed circular polarisation and \(r\) represents the mean reflectance under the filter for right handed circular polarisation.

Elytra

First export data frame with the Raw polarisation values in case it is needed

write.csv(RGB_ave, here::here("Data/FromCode/PolarzElytra_LERGBValues.csv"))

Thus we obtained the polarisation values for the elytron ROI by individual and by species:

# Calculate the polarisation:
ElyPol <-
  RGB_ave %>%
  spread(Type, RGB_ave) %>%
  dplyr::mutate(Polarization = (LCP - RCP) / (LCP + RCP))

head(ElyPol) # by individual

## # A tibble: 6 × 5
## # Groups:   ind [6]
##   ind      LCP   RCP   VIS Polarization
##   <chr>  <dbl> <dbl> <dbl>        <dbl>
## 1 abno01  7.40  6.61  5.31       0.0563
## 2 abno02  5.95  5.15  4.15       0.0723
## 3 anom01  4.18  3.17  2.52       0.136 
## 4 anom02  3.81  2.62  2.43       0.185 
## 5 atki01  7.53  4.00  5.17       0.306 
## 6 atki02  7.04  3.76  4.85       0.304

# To find the mean value per spp: 

ElyPolbb <-
  ElyPol %>%
  dplyr::select(ind, Polarization) %>%
  dplyr::mutate(spp = substr(ind, 1, 4)) # Add column for spp

ElyPolbb <-  as.data.frame(ElyPolbb)

# Import species names
SppNames <- read.csv(here::here("Data/9_CodesAndSpecies.csv"))

SppNames <-
  SppNames %>% 
  mutate (spp = substr(Ind, 1, 4)) %>% 
  filter (spp != "ambl"& # Amblyterus cicatricosus
          spp != "psqz"& # Pseudoschizongnatus schoenfeldti
          spp != "saul"& # Saulostomus villosus
          spp != "sqzb"& # Schizognathus burmeisteri
          spp != "sqzc"& # Schizognathus compressicornis
          spp != "sqzm"  # Schizognathus mesosternalis
            ) %>% # These species were removed. No phylogenetic info
  dplyr::select (-spp) %>% 
  dplyr::arrange(Ind) %>% 
  dplyr::rename("ind" = Ind)

# Combine the two data frames
ElyPolbyindN <- merge(ElyPolbb, SppNames, by="ind") # merge as "inner join"

# mean by species
EPolbs <-
  ElyPolbyindN%>%
  select(Polarization, phylogeny_name) %>% 
  dplyr::group_by(phylogeny_name) %>%
  dplyr::summarise(Pol = mean(Polarization)) # by species

Pronotum

Thus we obtained the polarisation values for the pronotum ROI by individual:

# Calculate the polarisation:
ProPol <-
  RGB_avePr %>%
  spread(Type, RGB_aveP) %>%
  dplyr::mutate(Polarization = (LCP - RCP) / (LCP + RCP))

head(ProPol) #  by individual

## # A tibble: 6 × 5
## # Groups:   ind [6]
##   ind      LCP   RCP   VIS Polarization
##   <chr>  <dbl> <dbl> <dbl>        <dbl>
## 1 abno01  6.60  5.68  4.68       0.0749
## 2 abno02  5.76  4.81  4.01       0.0891
## 3 anom01  4.15  3.01  2.60       0.160 
## 4 anom02  4.30  3.12  2.88       0.159 
## 5 atki01  3.28  2.65  2.41       0.107 
## 6 atki02  3.60  2.75  2.56       0.135

And by species:

# Average: One polarisation value by spp
ProPolbysp <-
  ProPol %>%
  dplyr::select(ind, Polarization) %>%
  dplyr::mutate(spp = substr(ind, 1, 4)) # Add column for spp

PPolbs <-
  ProPolbysp[2:3] %>%
  dplyr::group_by(spp) %>%
  dplyr::summarise(Pol = mean(Polarization)) # by species

head(PPolbs)

spp	Pol
abno	0.082
anom	0.16
atki	0.121
aurs	0.41
boid	0.116
brev	0.0606

Combined

Join the two data frames of polarisation by elytra and pronotum.

First make sure the two data frames contain the same number of species, thus the differences between the individuals in both data frames should be 0

ProPol <-
  ProPol %>%
  dplyr::rename(PronotPol = Polarization)

ElyPol <-
  ElyPol %>%
  dplyr::rename(ElytrPol = Polarization) # Keep all spp for other tests


ElyPolEP <-
  ElyPol %>%
  filter(ind != "ecry03" &
    ind != "mcla02" &
    ind != "roci02") # Removed, pronotum unavailable

inelytra <- levels(as.factor(ElyPolEP$ind))
inpronotum <- levels(as.factor(ProPol$ind))

setdiff(inelytra, inpronotum)

## character(0)

setdiff(inpronotum, inelytra)

## character(0)

Merge

# simplify Data frames
Elysimple <-
  ElyPolEP %>%
  select(ind, ElytrPol)

Prosimple <-
  ProPol %>%
  select(ind, PronotPol)

BeetlePol <-
  merge(Elysimple, Prosimple) %>% # Combine the two data frames
  dplyr::mutate(Difference = ElytrPol - PronotPol)

head(BeetlePol)

ind	ElytrPol	PronotPol	Difference
abno01	0.0563	0.0749	-0.0185
abno02	0.0723	0.0891	-0.0168
anom01	0.136	0.16	-0.0235
anom02	0.185	0.159	0.0262
atki01	0.306	0.107	0.199
atki02	0.304	0.135	0.169

We merged the polarisation data with the data frames obtained in previous steps that contain ecological variables, reflectivity, size and the species name as follows:

First compare the dimensions of the two data frames (Reflectance and polarisation):

Cons1 <-
  Cons1 %>%
  mutate(spp = substr(ind, 1, 4))

inTree <- levels(as.factor(Cons1$spp))
inRGB <- levels(as.factor(EPolbs$spp))

## Warning: Unknown or uninitialised column: `spp`.

setdiff(inTree, inRGB)

##  [1] "abno" "anom" "atki" "aurs" "boid" "brev" "brun" "clor" "conc" "ecry" "fchi" "gray" "hirs" "lats" "macl" "mcla" "meli" "mima" "mont" "mult" "narm" "neba" "nebn" "neus"
## [25] "oliv" "opal" "pali" "pczo" "pczp" "pczv" "pind" "pine" "poro" "prsi" "punc" "pvul" "rayn" "repa" "reps" "rept" "rina" "roci" "rose" "rubi" "rugo" "smgg" "smgp" "smgr"
## [49] "tars" "velu" "vrdi" "xyle" "xyls"

setdiff(inRGB, inTree)

## character(0)

Both contain the same species

# Compare length:
length(Cons1$ind)

## [1] 261

length(ElyPol$ind)

## [1] 177

# subset the consolidated by the ind in polariz data frame
Refbyind <- Cons1[Cons1$ind %in% ElyPol$ind, ]

But since they had different numbers of individuals we only considered the ones that were common for both data frames to combine them into one larger data frame

# Simplify the polarisation data
ElyPolbyind <-
  ElyPol %>%
  select(1, 5) # keep only one column and the individual

RPData <-
  merge(Refbyind, ElyPolbyind) # Combine the two

head(RPData)

ind	R_ALL	R_VIS	R_NIR	Res	size	phylogeny_name	PC1	PC2	spp	ElytrPol
abno01	24.4	15.2	32.4	-9.83	1.54	Anoplognathus_abnormis	-2.14	0.0716	abno	0.0563
abno02	30.2	17.9	40.9	-3.92	1.61	Anoplognathus_abnormis	-1.82	0.161	abno	0.0723
anom01	17.3	2.23	30.3	0.886	1.76	Anomala_antigua	-2.3	-2.28	anom	0.136
anom02	19.9	7.89	30.4	-4.6	1.77	Anomala_antigua	-0.386	-1.15	anom	0.185
atki01	15.5	6.19	23.6	-9.73	2.46	Calloodes_atkinsonii	-2.74	1.08	atki	0.306
atki02	18.7	6.64	29.2	-4.54	2.02	Calloodes_atkinsonii	-2.74	1.08	atki	0.304

We also summarised this data by species

RPDataAgg <-
  RPData %>%
  dplyr::select(-ind, -spp) %>% # remove individual id
  dplyr::select(phylogeny_name, everything()) %>% # order columns
  dplyr::group_by(phylogeny_name) %>% # group
  dplyr::summarise(across(everything(), list(mean))) # mean

head(RPDataAgg)

phylogeny_name	R_ALL_1	R_VIS_1	R_NIR_1	Res_1	size_1	PC1_1	PC2_1	ElytrPol_1
Anomala_antigua	18.6	5.06	30.4	-1.86	1.77	-1.34	-1.71	0.161
Anoplognathus_abnormis	27.3	16.5	36.6	-6.87	1.58	-1.98	0.116	0.0643
Anoplognathus_aeneus	22.1	10.4	32.3	-5.22	3.11	-1.9	0.419	0.112
Anoplognathus_aureus	42.6	34.7	49.5	-11.9	1.71	-3.49	-0.434	0.206
Anoplognathus_boisduvalii	31.1	17.3	43.2	-1.06	2.61	-1.64	-0.922	0.0885
Anoplognathus_brevicollis	27.4	12.9	40.1	0.236	2.51	-1.15	-3.63	0.0871

This way we can conduct the analysis at both individual and species level

Results from this step

As a general conclusion it does seem that there can be three mechanisms: White underlay, chiral and broadband chiral.

We obtained a data frame useful to compare the degree of polarisation between two body parts, elytron and pronotum

write.csv(BeetlePol, here::here("Data/FromCode/PolarzElytraPronot.csv"))

We also obtained data frames with the degree of polarisation by individual and by species:

write.csv(ElyPolbyindN, here::here("Data/FromCode/PolarzElytraByInd.csv")) # by individual
write.csv(EPolbs, here::here("Data/FromCode/PolarzElytraBySpp.csv")) # by species

Calculating degree of polarisation

Setting up

General

Parameters

Grey Standard

Raw RGB values

Calibrations

Elytra

Linearise

Equalise

Average

Pronotum

Linearise

Equalise

Average

Polarisation

Elytra

Pronotum

Combined

Results from this step