3d-shapes.qmd

---
jupyter: julia-1.10

execute:
  freeze: true  # re-render only when source changes
  cache: true
  warning: false
---

# 3d shape classification

In this lesson, we will explore a 3d-shape classification problem using some experimental approaches.

## The dataset

To get in touch with the classics, we are going to use the 3d shape dataset seen in [@singh2007topological].

![In [@singh2007topological] the authors used two different methods to cluster the different 3d shapes: an approximation to the Gromov Hausdorff distance, and a distance on the Mapper graph of each shape. Here in this chapter, we are going to use persistent homology.](images/3d-shapes/mapper.png)

The dataset can be found at [http://people.csail.mit.edu/sumner/research/deftransfer/data.html](http://people.csail.mit.edu/sumner/research/deftransfer/data.html). It consists of some reference pose of animals and human faces. The files used in this lesson can be found [in this Google drive](https://drive.google.com/drive/folders/1YCx80qTvyDvunPJ8f6Ae8w6tdRz0QLy2?usp=sharing). After downloading it, unzip the file and put then inside a directory called "meshes".

The files are written in the ".obj" format. They are *meshes*: sets of points and triangles that form a 3d image like the ones we can see in videogames.

Let's start with a flamingo shape. We load some libraries

```{julia}
# read meshes and plot
using Meshes, GeoIO
import CairoMakie as gl

# see progress
using ProgressMeter

# dataframes
using DataFramesMeta, CSV, Chain

# metric spaces and graphs
using MetricSpaces
using Graphs, SimpleWeightedGraphs

# persistent homology
import Ripserer
import PersistenceDiagrams as Pd
import Plots

# comparing the distance matrix
using Clustering, StatsPlots
import StatisticalMeasures.ConfusionMatrices as CM
using MultivariateStats
```

and define functions to read and visualize shapes

```{julia}
#| code-fold: true
function list_files(path="", pattern="")
    files =
        @chain begin
            map(walkdir(path)) do (root, dirs, files)
                joinpath.(root, files)
            end
            reduce(vcat, _)
            filter(x -> occursin(pattern, x), _)
        end

    files
end;

function get_class_files(class)
    @chain begin list_files("meshes/", class)
        filter(x -> occursin(".obj", x), _)
    end
end;

get_class_from_file(x) = split(x, "-")[1]

read_mesh(filepath) = GeoIO.load(filepath).geometry;

plot_mesh(ms) = viz(ms, showfacets = true, alpha = 0.5);
```

The reference pose is the following:

```{julia}
filepath = "meshes/flamingo-poses/flam-reference.obj"
ms = read_mesh(filepath)
plot_mesh(ms)
```

We can also visualize it as a gif

```{julia}
#| eval: false

# fig, ax, plt = plot_mesh(ms);

# gl.record(fig, "images/3d-shapes/flamingo-ref.gif", 0:31) do i
#     gl.rotate!(fig.scene, Vec3f(0, 0, 1), i * pi/16)
# end;
```

![](images/3d-shapes/flamingo-ref.gif)

and see its details:

![](images/3d-shapes/flamingo-details.png)


## Plotting {#sec-3d-plotting}

Here we will see all poses together, for each class. The last pose is the reference pose.

```{julia}
#| code-fold: true
function plot_entire_class(class) 
    files = get_class_files(class)
    fig = gl.Figure(size = (800, 1200));

    for (i, file) ∈ enumerate(files)
        p, q = divrem(i - 1, 3) .+ 1    
        ms = read_mesh(file)
        viz(fig[p, q], ms, title = "a")    
    end

    fig
end;

function create_gif_entire_class(class)
    fig = plot_entire_class(class)

    gl.record(fig, "images/3d-shapes/" * class * ".gif", 0:31) do i
        gl.rotate!(fig.scene, Vec3f(0, 0, 1), i * pi/16)
    end
end;
```

```{julia}
#| code-fold: true

# all_classes = [
#     "camel", "cat", "elephant", "face"
#     , "flamingo", "head", "horse", "lion"
#     ]

# @showprogress map(create_gif_entire_class, all_classes)
```

### Camel

![Camel poses](images/3d-shapes/camel.gif)

### Cat

![Cat poses](images/3d-shapes/cat.gif)

### Elephant

![Elephant poses](images/3d-shapes/elephant.gif)

### Face

![Face poses](images/3d-shapes/face.gif)

### Flamingo

![Flamingo poses](images/3d-shapes/flamingo.gif)

### Head

![Head poses](images/3d-shapes/head.gif)

### Horse

![Horse poses](images/3d-shapes/horse.gif)

### Lion

![Lion poses](images/3d-shapes/lion.gif)

## Setting the classification problem

We have 84 shapes in the following directories:

```{julia}
filter(!isfile, readdir("meshes/", join = true))
```

Each shape $s \in S$ has a class of the type camel, cat, elephant, etc. We can think of these classes as a function $c: S \to C$ where $C$ is the set of classes. Let $S_{rp}$ be the set of reference poses.

We will try to solve the following problem: can we correctly calculate $c(s)$ when we only know $c$ for $s \in S_{rp}$? That is: **knowing only the class of each reference pose, can we deduce the class of the remaining shapes using only the mesh file**?

In other words:

::: {.callout-note title="Problem"}

If I give you a 3d shape, can you say its class _if all you know is the reference pose of each class_?

:::

```{julia}
#| include: false

# function plot_reference_classes()
#     files = list_files("meshes/", "reference.obj")
#     fig = gl.Figure(size = (800, 1200));

#     for (i, file) ∈ enumerate(files)
#         p, q = divrem(i - 1, 3) .+ 1    
#         ms = read_mesh(file)
#         viz(fig[p, q], ms, title = "a")    
#     end

#     gl.record(fig, "images/3d-shapes/reference" * ".gif", 0:31) do i
#         gl.rotate!(fig.scene, Vec3f(0, 0, 1), i * pi/16)
#     end
# end
```

![The reference pose of each class.](images/3d-shapes/reference.gif)

This kind of problem is common in data science and is known as a "classification problem": we are trying to *atribute classes to objects, knowing the class of fewer other objects*.

## From meshes to metric spaces

::: {.callout-note title="An ingenuous attempt"}

An ingenuous attempt to solve the classification problem can be summarised as follows:

1. For each shape $S_i$, extract the points $X_i \subset \mathbb{R}^3$ and consider $d$ as the Euclidean distance;
1. Calculate the persistence diagram $D_i = dgm(X_i)$;
1. For each $D_i$, calculate the bottleneck distance from $D_i$ to all $D_j$ where $S_j$ is a reference pose;
1. The closest reference pose to $D_i$ will be the class of $S_i$.
:::

**This approach won't work** because of the two first steps:

- The euclidean distance is not appropriate for this problem. Flamingos in different poses will have a big Gromov-Hausdorff distance. We need to use some kind of geodesic distance.

- The amount of points in $X_i$ is too big to calculate the persistence diagram. The elephant class has more than 75.000 points for each shape. This will probably explode your RAM memory when calculating the Rips complex.

Fortunately, there are ways to contourn these problems!

- Extract a subset of "reasonably spaced points" of $S$ that still contains its core geometric properties;

- Calculate the geodesic distance between these points using the shape $S$.

### From meshes to $\mathbb{R}^3$

Let's extract the points of some $S$ as a subspace of $\mathbb{R}^3$:

```{julia}
mesh_to_metric_space(ms) = ms.vertices .|> coordinates .|> Vector |> EuclideanSpace;
```

```{julia}
X = mesh_to_metric_space(ms)
```

We can see that $X$ is made of 26907 points of $\mathbb{R}^3$. We can plot it:

```{julia}
gl.scatter(X, markersize = 1)
```

### From meshes to graphs

Now, to calculate the geodesic distance, we will create a graph from the mesh $S$. We define a function

```{julia}
#| code-fold: true
function graph_from_mesh(ms)
    # the set of vertices
    V = coordinates.(ms.vertices)

    # create an empty graph
    g = SimpleWeightedGraph()

    # add n vertices to it
    n = length(V)
    add_vertices!(g, n)

    # the set of triangles of the mesh ms
    triangles = ms.topology.connec

    # for each triangle, add its edges to the graph
    @showprogress desc = "Adding vertices to graph..." for t ∈ triangles
        v1, v2, v3 = t.indices
        add_edge!(g, v1, v2, dist_euclidean(V[v1], V[v2]))
        add_edge!(g, v1, v3, dist_euclidean(V[v1], V[v3]))
        add_edge!(g, v2, v3, dist_euclidean(V[v2], V[v3]))
    end

    g
end;
```

and create the graph $g$ from the mesh

```{julia}
g = graph_from_mesh(ms)
```

This weighted graph is the 1-skeleton of the mesh, and the weights between the vertices are the euclidean distance between then (as subsets of $\mathbb{R}^3$).

We can see the sparse array of its weight as follows:

```{julia}
weights(g)
```

Notice, however, that the mesh is not connected! This can be seen with

```{julia}
is_connected(g)
```

These are the connected components of $g$:

```{julia}
connected_components(g)
```

There is one big connected components, and several smaller ones with 1 point each. Let's extract the one with the most points and throw away the points of $X$ outside it.

```{julia}
#| code-fold: true

function extract_biggest_connected_component(g)
    cc_components = connected_components(g)
    ids_biggest_component = cc_components[findmax(length, cc_components)[2]]

    # modify the graph g on place
    g = g[ids_biggest_component]

    # return g and the ids of the biggest connected component
    g, ids_biggest_component
end;
```

```{julia}
g, ids_biggest_component = extract_biggest_connected_component(g);
```

We can see that $g$ now is connected:

```{julia}
is_connected(g)
```

Let's throw away from $X$ the points outside this component:

```{julia}
X = X[ids_biggest_component]
```

We now have 26394 points, which is a small reduction.

### Farthest points sampling

We could just select a random sample of points from our space, but points in high-density areas would be selected a lot more. There is another way to select points in a "well spaced manner", called the *farthest point sampling algorithm*. This algorithm was shown to me by Facundo Mémoli on a dirty blackboard in 2018 and the simplicity of it astonished me. For the curious ones, the algorithm is detailed below.

![An example of farthest point sample. The blue dots are the space $X$; the red dots are chosen each time far away from all the other red dots.](images/3d-shapes/farthest%20points.gif)

::: {.callout-important}
<!-- ## Algorithm -->

Let $(X, d)$ be a metric space. Fix an integer $n$. Let $C = \emptyset$ be the "set of chosen points". Select $x_1 \in X$ randomly and add it to $C$. Repeat the following procedure until you have $n$ points in $C$:

- Calculate the point $x \in X$ that is the most distant from all elements of $C$.

- Add $x$ to $C$.

- If $C$ has $n$ points, stop.

The set $C$ is called a *farthest points sampling* of $X$ with size $n$.

Notice that running the algorithm several times can lead to different sets $C$ because the first term is chosen randomly.
:::

Let's extract 400 points with the FPS algorithm and the euclidean distance:

```{julia}
ids_fps = farthest_points_sample(X, 400);
X_fps = X[ids_fps]
```


```{julia}
gl.scatter(X_fps, markersize = 10)
```

This is a very good approximation!

We are now interested in calculating the geodesic distance between these 400 points. But be careful! The geodesic distance need the entire mesh to work.

### Geodesic distances

Given a shape $S$, we can think of the geodesic distance between two points as "the least distance an ant would need to walk from one point to another". We will approximate this "walkable" paths using the edges of the triangles of the shape $S$. Remember: a mesh is a set of points and triangles!

The [Dijkstra algorithm](https://en.wikipedia.org/wiki/Dijkstra%27s_algorithm) is perfect for our needs: it calculates the shortest path from one point to another in a weighted graph. So all we need is to:

- Transform $S$ into a graph where the edges have weights (the euclidean distance between these points);
- Calculate the shortest path between each two points.

We already have the first item, so let's calculate the second.

```{julia}
#| code-fold: true

function geodesic_distance_from_graph(g, ids)
    n = length(ids)
    D = zeros(n, n)

    # for each point, calculate the distance from it to every other point of g
    @showprogress desc = "Calculating geodesic distance..." Threads.@threads for (i, id) ∈ collect(enumerate(ids))
        dts = dijkstra_shortest_paths(g, id)
        D[i, :] = dts.dists[ids]
    end

    # force simmetry on X, because of small difference
    # in the calculation of paths
    for i ∈ 1:n
        for j ∈ i:n
            D[i, j] = D[j, i]
        end
    end

    # normalize the distance so the max is 1
    max_dist = maximum(D)
    D = D ./ max_dist

    return D
end;
```

```{julia}
D = geodesic_distance_from_graph(g, ids_fps)
```

We can see that $D$ makes sense just by plotting $X_fps$ colored by the sum of the distances to each points:

```{julia}
exc = map(sum, eachcol(D))

gl.scatter(X_fps, color = exc, markersize = 10)
```

Looks good! The extremities of the flamingo are in a lighter color, indicating that the sum of the distances there is bigger. Now we have 1000 points sampled from $S$, together with the geodesic distance.

## Persistent homology

We can now calculate the persistence diagram of $X_{fps}$ with the geodesic distance and use it! Let's load some packages and calculate it

```{julia}
pd = Ripserer.ripserer(D, dim_max = 2, verbose=true, sparse = true, threshold = 0.4)
```

Ploting the intervals looks as follows:

```{julia}
#| code-fold: true

function plot_barcode(pd)
    # get the size of the longest interval
    threshold = 
        @chain begin
            vcat(pd...)
            last.(_)
            filter(isfinite, _)
            maximum
        end

    # plot the barcode using this interval as the maximum value of the x-axis
    Ripserer.barcode(pd, infinity = threshold)
end;
```

```{julia}
plot_barcode(pd)
```

or just the 1- and 2-dimensional barcode:

```{julia}
plot_barcode(pd[2:3])
```

## Summarizing

All the hard work on the previous sections was just to prepare our dataset from file to barcode. That's why they say that **data science is 80% preparing the data and 20% analyzing it!**

We can summarise what we did with the following function:


```{julia}
#| code-fold: true

function file_to_barcode(filepath; n_points = 1000, dim_max = 1)
    ms = read_mesh(filepath)

    X = mesh_to_metric_space(ms)
    g = graph_from_mesh(ms)

    g, ids_biggest_component = extract_biggest_connected_component(g)
    X = X[ids_biggest_component]

    ids_fps = farthest_points_sample(X, n_points);
    X_fps = X[ids_fps]

    D = geodesic_distance_from_graph(g, ids_fps)

    pd = Ripserer.ripserer(D, dim_max = dim_max, verbose=true, sparse = true, threshold = 0.8)

    return X_fps, D, pd
end;
```

We also define some functions to save the barcodes and metric spaces to disk, so we don't have to calculate all of them in a single session. Calculating the 2-dimensional barcode can take some time depending on your hardware!

```{julia}
#| code-fold: true
function pd_to_dataframe(pd)
    df = @chain begin
        map(pd) do p
            DataFrame(
                birth=p .|> first, death=p .|> last, dim=p.dim
            )
        end
        vcat(_...)
    end

    df
end;

function dataframe_to_pd(df)
    df.threshold .= 1

    @chain df begin
        groupby(:dim)
        collect
        map(Pd.PersistenceDiagram, _)
    end
end;

function metric_space_to_df(X) 
    @chain X_fps begin
        stack
        transpose
        DataFrame(_, :auto)
    end
end;
```

Now we loop over all meshes, calculate its persistence diagram and save it to disk, together with the $X_{fps}$ metric space as above.

**Important**: This can take some time! If you cloned my github repository, these files are already there, so you can skip the following piece of code:

```{julia}
#| eval: false
overwrite_old_files = true

@showprogress "Calculating barcode..." for file ∈ list_files("meshes/", ".obj") 
    println("Calculating barcode from file $file ...")

    file_pd = replace(file, ".obj" => "-pd.csv")
    # skip if there is a file already
    if isfile(file_pd) & !overwrite_old_files continue end

    X_fps, D, pd = file_to_barcode(file, n_points = 350, dim_max = 2)
    df = pd_to_dataframe(pd)

    CSV.write(file_pd, df)

    file_X = replace(file, ".obj" => "-points.csv")
    CSV.write(file_X, metric_space_to_df(X_fps))
end
```

We read the persistence diagrams saved on disk and pass them to table (a DataFrame object), but first we throw away small intervals.

```{julia}
#| code-fold: true

function throw_away_small_intervals(pd, min_pers = 0.01)
    map(pd) do p
        filter(x -> Pd.persistence(x) > min_pers, p)
    end
end;

function read_pds_from_files(directory, min_interval_size = 0.05)
  pds_df = DataFrame()

#   file = list_files("meshes/", "-pd.csv")[1]
  for file ∈ list_files(directory, "-pd.csv")

      pd = @chain begin
          CSV.read(file, DataFrame)
          dataframe_to_pd(_)
          throw_away_small_intervals(min_interval_size)
      end

      name = replace(file, "-pd.csv" => "")
      push!(pds_df, (Path = name, Persistence_diagram = pd))
  end

  pds_df

  sort!(pds_df, :Path)

  pds_df.File = [split(s, "/")[3] for s ∈ pds_df.Path]
  pds_df.Class = [split(s, "-")[1] for s ∈ pds_df.File]

  pds_df
end;
```

## A metric problem

Now we proceed to create a dataframe that contains all we need to classify our shapes.

```{julia}
pds_df = read_pds_from_files("meshes/", 0.01);

first(pds_df, 20)
```

The dataframe looks ok! You can plot the barcodes as follows:

```{julia}
pd2 = pds_df.Persistence_diagram[1]
plot_barcode(pd2[2:3])
```

Now we calculate the bootleneck distance between each pair of persistence diagrams. This can take some time! If you cloned the repository, you don't need to run this piece of code.

```{julia}
#| eval: false
pds = pds_df.Persistence_diagram

n = nrow(pds_df)
DB = zeros(n, n)

@showprogress for i ∈ 1:n
    for j ∈ i:n
        if i == j
            DB[i, j] = 0 
            continue 
        end

        DB[i, j] = 
            Pd.Bottleneck()(pds[i][2], pds[j][2]) + 
            Pd.Bottleneck()(pds[i][3], pds[j][3])

        DB[j, i] = DB[i, j]
    end
end

CSV.write("meshes/bottleneck_distance.csv", DataFrame(DB, :auto))
```

Notice that we defined the distance $DB_{i, j}$ between two shapes $X_i$ and $X_j$ as 

$$
DB_{i, j} = d_b(\text{dgm}_1(X_i), \text{dgm}_1(X_j)) + d_b(\text{dgm}_2(X_i), \text{dgm}_2(X_j))
$$

where $d_b$ is the bottleneck distance, and $dgm_i$ is the $i$-dimensional persistence diagram. Notice that **we did not use the 0-dimensional persistent homology, because all shapes are connected**.

## Facing the truth

We read $DB$ from disk, in case you did not calculate it previously

```{julia}
DB = CSV.read("meshes/bottleneck_distance.csv", DataFrame) |> Matrix
DB
```

How can we visualize $DB$?

### Heatmap

Heatmaps are matrices (2d arrays) where each element is colored by some value. In the ideal world, our matrix $DB$ will have a square of low values for each comparison between the same class, and a higher value when comparing elements of different classes.

```{julia}
labels = pds_df.Class |> copy
for i ∈ 2:length(pds_df.Class)
    if labels[i] == pds_df.Class[i-1]
        labels[i] = ""
    end
end

labels = (1:nrow(pds_df), labels)

plot(DB, st = :heatmap, xticks = labels, yticks = labels)
```

We can see that the diagonals (intra-class comparisons) are good enough: they are dark-coloured squares, indicating small distances. But there are also black squares outside the diagonal, which means that some different classes are close to each other. 

### Dendrograms

Dendrograms are a nice way to visualize a distance matrix. They represent the evolution of the connected components as some parameter (usually the radius of circles centered on each point) grows.

```{julia}
#| code-fold: true
function plot_hc(hc)
    plot(
        hc, xticks = (1:nrow(pds_df), pds_df.File[hc.order])
        , xflip = true, xrotation = 270
        , xtickfont = font(5, "Roboto")
        )
end;
```

The [single linkage dendrogram](https://en.wikipedia.org/wiki/Single-linkage_clustering) is a mess:

```{julia}
hclust(DB, linkage=:single, branchorder = :optimal) |> plot_hc
```

But the [complete linkage dendrogram](https://en.wikipedia.org/wiki/Complete-linkage_clustering) shows more hope:

```{julia}
hclust(DB, linkage = :complete, branchorder = :optimal) |> plot_hc
```

Here is also the Ward algorithm:

```{julia}
hclust(DB, linkage = :ward, branchorder = :optimal) |> plot_hc
```

### Accuracy

Now let's calculate for each shape which reference class is closer to it:

```{julia}
#| code-fold: true

function get_score_dataframe(pds_df, DB) 
  score = @select(pds_df, :Path, :File, :Class)

  score.Nearest_class .= ""

  ids_reference = findall(x -> occursin("-reference", x), score.Path)
  names_reference = score.File[ids_reference]

  for i ∈ 1:nrow(pds_df)
      id = sortperm(DB[i, ids_reference])[1]
      score.Nearest_class[i] = names_reference[id]
  end

  score.Right_class =
      get_class_from_file.(score.Class) .== 
      get_class_from_file.(score.Nearest_class)

  score
end;
```

```{julia}
score = get_score_dataframe(pds_df, DB)
```

Our accuracy was

```{julia}
score.Right_class |> mean
```

that is: 82%. Not bad, but not excellent either!

Here is the [confusion matrix](https://en.wikipedia.org/wiki/Confusion_matrix):

```{julia}
score_new = @rsubset score !occursin("ref", :File)
x = CM.confmat(get_class_from_file.(score_new.Nearest_class), score_new.Class)
```

![Confusion matrix for DB.](images/3d-shapes/confusion%20matrix.png)

Some conclusions we can take from the above analysis:

- Camels, elephants and faces are all correct;
- There is one head mistaken as a face;
- There are two mistaken horses;
- Cats are confused with lions.

Why did this happen? Look back at the rotating plots in section [@sec-3d-plotting]: a cat and a lion are almost isometric in the geodesic sense! Thus, **no tool that uses only the metric will be able to split these classes**.

### MDS plot

Another last visualization technique is the [multidimensional scaling](https://en.wikipedia.org/wiki/Multidimensional_scaling) plot. This method takes a distance matrix and tries to project it down to $\mathbb{R}^2$) (or $\mathbb{R}^3$) while distorting the distances as little as possible.

```{julia}
#| code-fold: true

function mds_plot(D, score)
    M = fit(MDS, D; distances = true, maxoutdim = 2)
    Y = predict(M)

    score.Row = 1:length(score.Class)

    dfs = @chain score begin
        groupby(:Class)
        collect
    end

    fig = gl.Figure();
    ax = gl.Axis(fig[1,1])

    colors = cgrad(:tableau_10, 8, categorical = true)

    for (i, df) ∈ enumerate(dfs)    
        gl.scatter!(
            ax, Y[:, df.Row]
            , label = df.Class[1], markersize = 15
            , color = colors[i]
            )
    end

    gl.axislegend();
    fig

    fig
end
```

```{julia}
mds_plot(DB, score)
```

Can we do better?

With only the distance matrix $DB$, as we could see above, the answer is **NO**. We need to find more tools to compare these shapes.

## Using volumes

::: {.callout-caution title="Warning"}

Don't try the following scenes at home! Read the text and explain why this approach is wrong.

:::


Each mesh has a notion of volume that can be calculated with the function `measure`.

We can see below the mean volume for each class:

```{julia}
volumes = @showprogress map(score.Path) do f
    ms = read_mesh(f * ".obj")
    measure(ms)
end

df_volumes = DataFrame(
    Path = score.Path, Class = score.Class, Volume = volumes
    )

@chain df_volumes begin
    @groupby :Class
    @combine :Mean_volume = mean(:Volume)
end
```

The volume is enough to separate all the classes except lions and flamingos (which have volumes close to 0.54).

So let's add this difference in volume to our original matrix $DB$, and call it $DB2$:

```{julia}
n = nrow(pds_df)
volume_difs = zeros(n, n)

for i ∈ 1:n
    for j ∈ i:n
        volume_difs[i, j] = abs(df_volumes.Volume[i] - df_volumes.Volume[j])
        volume_difs[j, i] = volume_difs[i, j]
    end
end

# ignore differences in volume greater than maximum(DB)
replace!(x -> x > maximum(DB) ? maximum(DB) : x, volume_difs)

DB2 = DB .+ volume_difs
```

Its dendrogram is the following:

```{julia}
hclust(DB2, linkage = :complete, branchorder = :optimal) |> plot_hc
```

which is much better than before!

If again we calculate the closest reference pose to each given shape, we now get

```{julia}
score2 = get_score_dataframe(pds_df, DB2)
```

```{julia}
score2.Right_class |> mean
```

nearly 100%!

The only misclassification can be seen with

```{julia}
@rsubset score2 !:Right_class 
```

where a head was mistaken for a face.