-
Notifications
You must be signed in to change notification settings - Fork 0
/
README.Rmd
391 lines (314 loc) · 28.1 KB
/
README.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
---
output:
github_document:
pandoc_args: --webtex=https://ibm.codecogs.com/png.latex?
always_allow_html: yes
---
<!-- README.md is generated from README.Rmd. Please edit that file -->
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.path = "man/figures/README-",
out.width = "100%"
)
```
# cmcR
<!-- badges: start -->
[](https://app.codecov.io/gh/CSAFE-ISU/cmcR?branch=master/)
[](https://github.com/CSAFE-ISU/cmcR/actions/)
<!-- badges: end -->
The cmcR package provides an open-source implementation of the Congruent Matching Cells method for cartridge case identification as proposed by [Song (2013)](https://tsapps.nist.gov/publication/get_pdf.cfm?pub_id=911193) as well as the "High CMC" method proposed by [Tong et al. (2015)](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4730689/pdf/jres.120.008.pdf).
## Installation
<!-- You can install the released version of cmcR from [CRAN](https://CRAN.R-project.org) with: -->
<!-- ``` r -->
<!-- install.packages("cmcR") -->
<!-- ``` -->
Install the development version from [GitHub](https://github.com/CSAFE-ISU/cmcR) with:
``` r
# install.packages("devtools")
devtools::install_github("jzemmels/cmcR")
```
Cartridge case scan data can be accessed at the [NIST Ballisitics Toolmark Research Database](https://tsapps.nist.gov/NRBTD/Studies/Search)
## Example
We will illustrate the package's functionality here. Please refer to the package vignettes available under the "Articles" tab of the [package website](https://csafe-isu.github.io/cmcR/index.html) for more information.
```{r setup,warning=FALSE,message=FALSE}
library(cmcR)
library(magrittr)
library(dplyr)
library(ggplot2)
library(x3ptools)
```
Consider the known match cartridge case pair Fadul 1-1 and Fadul 1-2. The `read_x3p` function from the [x3ptools](https://github.com/heike/x3ptools) package can read scans from the [NBTRD](https://tsapps.nist.gov/NRBTD/Studies/Search) given the appropriate address. The two scans are read below and visualized using the [`x3pListPlot`](https://csafe-isu.github.io/cmcR/reference/x3pListPlot.html) function.
```{r}
fadul1.1_id <- "DownloadMeasurement/2d9cc51f-6f66-40a0-973a-a9292dbee36d"
# Same source comparison
fadul1.2_id <- "DownloadMeasurement/cb296c98-39f5-46eb-abff-320a2f5568e8"
# Code to download breech face impressions:
nbtrd_url <- "https://tsapps.nist.gov/NRBTD/Studies/CartridgeMeasurement/"
fadul1.1_raw <- x3p_read(paste0(nbtrd_url,fadul1.1_id)) %>%
x3ptools::x3p_scale_unit(scale_by = 1e6)
fadul1.2_raw <- x3p_read(paste0(nbtrd_url,fadul1.2_id)) %>%
x3ptools::x3p_scale_unit(scale_by = 1e6)
x3pListPlot(list("Fadul 1-1" = fadul1.1_raw,
"Fadul 1-2" = fadul1.2_raw),
type = "faceted")
```
### Preprocessing
To perform a proper comparison of these two cartridge cases, we need to remove regions that do not come into uniform or consistent contact with the breech face of the firearm. These include the small clusters of pixels in the corners of the two scans from the microscope staging area, and the plateaued region of points around the firing pin impression hole near the center of the scan. A variety of processing procedures are implemented in the cmcR package. Functions of the form `preProcess_*` perform the preprocessing procedures. See the [funtion reference](https://csafe-isu.github.io/cmcR/reference/index.html) of the cmcR package for more information regarding these procedures. As is commonly done when comparing cartridge cases, we downsample each scan (by a factor of 4, selecting every other row/column) using the `sample_x3p` function.
```{r}
fadul1.1_processed <- fadul1.1_raw %>%
preProcess_crop(region = "exterior",offset = -30) %>%
preProcess_crop(region = "interior",offset = 200) %>%
preProcess_removeTrend(statistic = "quantile",
tau = .5,
method = "fn") %>%
preProcess_gaussFilter() %>%
x3ptools::x3p_sample()
fadul1.2_processed <- fadul1.2_raw %>%
preProcess_crop(region = "exterior",offset = -30) %>%
preProcess_crop(region = "interior",offset = 200) %>%
preProcess_removeTrend(statistic = "quantile",
tau = .5,
method = "fn") %>%
preProcess_gaussFilter() %>%
x3ptools::x3p_sample()
x3pListPlot(list("Processed Fadul 1-1" = fadul1.1_processed,
"Processed Fadul1-2" = fadul1.2_processed),
type = "faceted")
```
### Cell-based comparison procedure
Functions of the form `comparison_*` perform the steps of the cell-based comparison procedure. The data generated from the cell-based comparison procedure are kept in a [`tibble`](https://tibble.tidyverse.org/) where one row represents a single cell/region pairing.
The `comparison_cellDivision` function divides a scan up into a grid of cells. The `cellIndex` column represents the `row,col` location in the original scan each cell inhabits. Each cell is stored as an `.x3p` object in the `cellHeightValues` column. The benefit of using a `tibble` structure is that processes such as removing rows can be accomplished using simple `dplyr` commands such as `filter`.
```{r}
cellTibble <- fadul1.1_processed %>%
comparison_cellDivision(numCells = c(8,8))
cellTibble
```
The `comparison_getTargetRegions` function extracts a region from a target scan (in this case Fadul 1-2) to be paired with each cell in the reference scan.
```{r}
cellTibble <- cellTibble %>%
mutate(regionHeightValues = comparison_getTargetRegions(cellHeightValues = cellHeightValues,
target = fadul1.2_processed))
cellTibble
```
We want to exclude cells and regions that are mostly missing from the scan. The `comparison_calcPropMissing` function calculates the proportion of missing values in a surface matrix. The call below excludes rows in which either the cell or region contain more that 85% missing values.
```{r}
cellTibble <- cellTibble %>%
mutate(cellPropMissing = comparison_calcPropMissing(cellHeightValues),
regionPropMissing = comparison_calcPropMissing(regionHeightValues)) %>%
filter(cellPropMissing <= .85 & regionPropMissing <= .85)
cellTibble %>%
select(cellIndex,cellPropMissing,regionPropMissing)
```
We can standardize the surface matrix height values by centering/scaling by desired functions (e.g., mean and standard deviation). Also, to apply frequency-domain techniques in comparing each cell and region, the missing values in each scan need to be replaced. These operations are performed in the `comparison_standardizeHeightValues` and `comparison_replaceMissingValues` functions.
Then, the `comparison_fft_ccf` function estimates the translations required to align the cell and region using the [Cross-Correlation Theorem](https://mathworld.wolfram.com/Cross-CorrelationTheorem.html). The `comparison_fft_ccf` function returns a data frame of 3 `x`, `y`, and `fft_ccf` values: the $x,y$ estimated translation values at which the CCF$_\max$ value is attained between the cell and region. The `tidyr::unnest` function can unpack the data frame into 3 separate columns, if desired.
```{r}
cellTibble <- cellTibble %>%
mutate(cellHeightValues = comparison_standardizeHeights(cellHeightValues),
regionHeightValues = comparison_standardizeHeights(regionHeightValues)) %>%
mutate(cellHeightValues_replaced = comparison_replaceMissing(cellHeightValues),
regionHeightValues_replaced = comparison_replaceMissing(regionHeightValues)) %>%
mutate(fft_ccf_df = comparison_fft_ccf(cellHeightValues = cellHeightValues_replaced,
regionHeightValues = regionHeightValues_replaced))
cellTibble %>%
tidyr::unnest(cols = fft_ccf_df) %>%
select(cellIndex,fft_ccf,x,y)
```
Because so many missing values need to be replaced, the CCF$_{\max}$ value calculated in the `fft_ccf` column using frequency-domain techniques is not a very good similarity score (doesn't differentiate matches from non-matches well).
However, the `x` and `y` estimated translations are good estimates of the "true" translation values needed to align the cell and region.
To calculate a more accurate similarity score, we can use the pairwise-complete correlation in which only pairs of non-missing pixels are considered in the correlation calculation.
To calculate this the `comparison_alignedTargetCell` function takes the cell, region, and CCF-based alignment information and returns a matrix of the same dimension as the reference cell representing the sub-matrix of the target region that the cell aligned to.
We can then calculate the pairwise-complete correlation as shown below.
```{r}
cellTibble %>%
dplyr::mutate(alignedTargetCell = comparison_alignedTargetCell(cellHeightValues = .data$cellHeightValues,
regionHeightValues = .data$regionHeightValues,
target = fadul1.2_processed,
theta = 0,
fft_ccf_df = .data$fft_ccf_df)) %>%
dplyr::mutate(pairwiseCompCor = purrr::map2_dbl(.data$cellHeightValues,.data$alignedTargetCell,
~ cor(c(.x$surface.matrix),c(.y$surface.matrix),
use = "pairwise.complete.obs"))) %>%
tidyr::unnest("fft_ccf_df") %>%
select(cellIndex,x,y,pairwiseCompCor)
```
Finally, this entire comparison procedure is to be repeated over a number of rotations of the target scan. The entire cell-based comparison procedure is wrapped in the `comparison_allTogether` function. The resulting data frame below contains the features that are used in the decision-rule procedure
```{r}
kmComparisonFeatures <- purrr::map_dfr(seq(-30,30,by = 3),
~ comparison_allTogether(reference = fadul1.1_processed,
target = fadul1.2_processed,
theta = .,
returnX3Ps = TRUE))
kmComparisonFeatures %>%
arrange(theta,cellIndex)
```
### Decision rule
The decision rules described in [Song (2013)](https://tsapps.nist.gov/publication/get_pdf.cfm?pub_id=911193) and [Tong et al. (2015)](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4730689/pdf/jres.120.008.pdf) are implemented via the `decision_*` functions.
We refer to the two decision rules as the original method of [Song (2013)](https://tsapps.nist.gov/publication/get_pdf.cfm?pub_id=911193) and the High CMC method, respectively.
Considering the `kmComparisonFeatures` data frame returned above, we can interpret both of these decision rules as logic that separates "aberrant" from "homogeneous" similarity features.
The two decision rules principally differ in how they define an homogeneity.
The original method of [Song (2013)](https://tsapps.nist.gov/publication/get_pdf.cfm?pub_id=911193) considers only the similarity features at which the maximum correlation is attained for each cell across all rotations considered. Since there is ambiguity in exactly how the correlation is computed in the original methods, we will consider the features at which specifically the maximum `pairwiseCompCor` is attained (instead of using the `fft_ccf` column).
```{r}
kmComparisonFeatures %>%
group_by(cellIndex) %>%
top_n(n = 1,wt = pairwiseCompCor)
```
The above set of features can be thought of as the `x`, `y`, and `theta` "votes" that each cell most strongly "believes" to be the correct alignment of the entire scan. If a pair is truly matching, we would expect many of these votes to be similar to each other; indicating that there is an approximate consensus of the true alignment of the entire scan (at least, this is the assumption made in [Song (2013)](https://tsapps.nist.gov/publication/get_pdf.cfm?pub_id=911193)). In [Song (2013)](https://tsapps.nist.gov/publication/get_pdf.cfm?pub_id=911193), the consensus is defined to be the median of the `x`, `y`, and `theta` values in this `topVotesPerCell` data frame. Cells that are deemed "close" to these consensus values and that have a "large" correlation value are declared Congruent Matching Cells (CMCs). Cells with `x`, `y`, and `theta` values that are within user-defined $T_x$, $T_y$, and $T_{\theta}$ thresholds of the consensus `x`, `y`, and `theta` values are considered "close." If these cells also have a correlation greater than a user-defined $T_{\text{CCF}}$ threshold, then they are considered CMCs. Note that these thresholds are chosen entirely by experimentation in the CMC literature.
```{r}
kmComparison_originalCMCs <- kmComparisonFeatures %>%
mutate(originalMethodClassif = decision_CMC(cellIndex = cellIndex,
x = x,
y = y,
theta = theta,
corr = pairwiseCompCor,
xThresh = 20,
yThresh = 20,
thetaThresh = 6,
corrThresh = .5))
kmComparison_originalCMCs %>%
filter(originalMethodClassif == "CMC")
```
Many have pointed out that there tends to be many cell/region pairs that exhibit high correlation other than at the "true" rotation (theta) value.
In particular, a cell/region pair may attain a very high correlation at the "true" theta value, yet attain its maximum correlation at a theta value far from the consensus theta value.
The original method of [Song (2013)](https://tsapps.nist.gov/publication/get_pdf.cfm?pub_id=911193) is quite sensitive to this behavior.
The original method of [Song (2013)](https://tsapps.nist.gov/publication/get_pdf.cfm?pub_id=911193) only considers the "top" vote of each cell/region pairing, so it is not sensitive to how that cell/region pairing behaves across multiple rotations.
[Tong et al. (2015)](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4730689/pdf/jres.120.008.pdf) propose a different decision rule procedure that considers the behavior of cell/region pairings across multiple rotations. This method would come to be called the High CMC method. The procedure involves computing a "CMC-`theta`" distribution where for each value of `theta`, the `x` and `y` values are compared to consensus `x` and `y` values (again, the median) and the correlation values to a minimum threshold. A cell is considered a "CMC candidate" (our language, not theirs) at a particular `theta` value if its `x` and `y` values are within $T_x$, $T_y$ thresholds of the consensus `x` and `y` values and the correlation is at least as larges as the $T_{\text{CCF}}$ threshold. This is similar to the original method of [Song (2013)](https://tsapps.nist.gov/publication/get_pdf.cfm?pub_id=911193) except that it relaxes the requirement that the top `theta` value be close to a consensus. Continuing with the voting analogy, think of this as an [approval voting system](https://en.wikipedia.org/wiki/Approval_voting) where each cell is allowed to vote for multiple `theta` values as long as the `x` and `y` votes are deemed close to the `theta`-specific `x`,`y` consensuses and the correlation values are sufficiently high.
The CMC-`theta` distribution consists of the "CMC candidates" at each value of `theta`. The assumption made in [Tong et al. (2015)](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4730689/pdf/jres.120.008.pdf) is that, for a truly matching cartridge case pair, a large number of CMC candidates should be concentrated around true `theta` alignment value. In their words, the CMC-`theta` distribution should exhibit a "prominent peak" close to the rotation at which the two cartridge cases actually align. Such a prominent peak should not occur for a non-match cartridge case pair.
The figure below shows an example of a CMC-`theta` distribution between Fadul 1-1 and Fadul 1-2 constructed using the `decision_highCMC_cmcThetaDistrib` function. We can clearly see that a mode is attained around -24 degrees.
```{r}
kmComparisonFeatures %>%
mutate(cmcThetaDistribClassif = decision_highCMC_cmcThetaDistrib(cellIndex = cellIndex,
x = x,
y = y,
theta = theta,
corr = pairwiseCompCor,
xThresh = 20,
yThresh = 20,
corrThresh = .5)) %>%
filter(cmcThetaDistribClassif == "CMC Candidate") %>%
ggplot(aes(x = theta)) +
geom_bar(stat = "count",
alpha = .7) +
theme_bw() +
ylab("CMC Candidate Count") +
xlab(expression(theta))
```
The next step of the High CMC method is to automatically determine if a mode (i.e., a "prominent peak") exists in a CMC-`theta` distribution. If we find a mode, then there is evidence that a "true" rotation exists to align the two cartridge cases implying the cartridge cases must be matches (such is the logic employed in [Tong et al. (2015)](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4730689/pdf/jres.120.008.pdf)). To automatically identify a mode, [Tong et al. (2015)](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4730689/pdf/jres.120.008.pdf) propose determining the range of `theta` values with "high" CMC candidate counts (if this range is small, then there is likely a mode). They define a "high" CMC candidate count to be $CMC_{\text{high}} \equiv CMC_{\max} - \tau$ where $CMC_{\max}$ is the maximum value attained in the CMC-`theta` distribution (17 in the plot shown above) and $\tau$ is a user-defined constant ([Tong et al. (2015)](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4730689/pdf/jres.120.008.pdf) use $\tau = 1$). Any `theta` value with associated an associated CMC candidate count at least as large as $CMC_{\text{high}}$ have a "high" CMC candidate count while any others have a "low" CMC candidate count.
The figure below shows the classification of `theta` values into "High" and "Low" CMC candidate count groups using the `decision_highCMC_identifyHighCMCThetas` function. The High CMC count threshold is shown as a dashed line at $21 - 1$ CMCs.
```{r}
kmComparisonFeatures %>%
mutate(cmcThetaDistribClassif = decision_highCMC_cmcThetaDistrib(cellIndex = cellIndex,
x = x,
y = y,
theta = theta,
corr = pairwiseCompCor,
xThresh = 20,
yThresh = 20,
corrThresh = .5)) %>%
decision_highCMC_identifyHighCMCThetas(tau = 1) %>%
filter(cmcThetaDistribClassif == "CMC Candidate") %>%
ggplot() +
geom_bar(aes(x = theta, fill = thetaCMCIdentif),
stat = "count",
alpha = .7) +
geom_hline(aes(yintercept = max(cmcCandidateCount) - 1),
linetype = "dashed") +
scale_fill_manual(values = c("black","gray50")) +
theme_bw() +
ylab("CMC Candidate Count") +
xlab(expression(theta))
```
If the range of High CMC count `theta` values is less than the user-defined $T_{\theta}$ threshold, then [Tong et al. (2015)](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4730689/pdf/jres.120.008.pdf) classify all CMC candidates in the identified `theta` mode as actual CMCs.
The `decision_CMC` function classifies CMCs based on this High CMC criterion if a value for `tau` is given. Note that it internally calls the `decision_highCMC_cmcThetaDistrib` and `decision_highCMC_identifyHighCMCThetas` functions (although they are exported as diagnostic tools). A cell may be counted as a CMC for multiple `theta` values. In these cases, we will only consider the alignment values at which the cell attained its maximum CCF and was classified as a CMC. If the cartridge case pair "fails" the High CMC criterion (i.e., the range of High CMC candidate `theta` values is deemed too large), every cell will be classified as "non-CMC (failed)" under the High CMC method. When it comes to combining the CMCs from two comparison directions (cartridge case A vs. B and B vs. A), we must treat a cell classified as a non-CMC because the High CMC criterion failed differently from a cell classified as a non-CMC for which the High CMC criterion passed.
<!-- [Tong et al. (2015)](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4730689/pdf/jres.120.008.pdf) propose using the CMC count determined under the original method of [Song (2013)](https://tsapps.nist.gov/publication/get_pdf.cfm?pub_id=911193) as a backup CMC count. -->
```{r}
kmComparison_highCMCs <- kmComparisonFeatures %>%
mutate(highCMCClassif = decision_CMC(cellIndex = cellIndex,
x = x,
y = y,
theta = theta,
corr = pairwiseCompCor,
xThresh = 20,
yThresh = 20,
thetaThresh = 6,
corrThresh = .5,
tau = 1))
#Example of cells classified as CMCs and non-CMCs
kmComparison_highCMCs %>%
slice(21:35)
```
In summary: the `decison_CMC` function applies either the decision rules of the original method of [Song (2013)](https://tsapps.nist.gov/publication/get_pdf.cfm?pub_id=911193) or the High CMC method of [Tong et al. (2015)](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4730689/pdf/jres.120.008.pdf), depending on whether the user specifies a value for the High CMC threshold `tau`.
```{r}
kmComparison_allCMCs <- kmComparisonFeatures %>%
mutate(originalMethodClassif = decision_CMC(cellIndex = cellIndex,
x = x,
y = y,
theta = theta,
corr = pairwiseCompCor,
xThresh = 20,
thetaThresh = 6,
corrThresh = .5),
highCMCClassif = decision_CMC(cellIndex = cellIndex,
x = x,
y = y,
theta = theta,
corr = pairwiseCompCor,
xThresh = 20,
thetaThresh = 6,
corrThresh = .5,
tau = 1))
#Example of cells classified as CMC under 1 decision rule but not the other.
kmComparison_allCMCs %>%
slice(21:35)
```
The set of CMCs computed above are based on assuming Fadul 1-1 as the reference scan and Fadul 1-2 as the target scan. [Tong et al. (2015)](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4730689/pdf/jres.120.008.pdf) propose performing the cell-based comparison and decision rule procedures with the roles reversed and combining the results. They indicate that if the High CMC method fails to identify a `theta` mode in the CMC-`theta` distribution, then the minimum of the two CMC counts computed under the original method of [Song (2013)](https://tsapps.nist.gov/publication/get_pdf.cfm?pub_id=911193) should be used as the CMC count (although they don't detail how to proceed if a `theta` mode is identified in only one of the two comparisons).
```{r}
#Compare using Fadul 1-2 as reference and Fadul 1-1 as target
kmComparisonFeatures_rev <- purrr::map_dfr(seq(-30,30,by = 3),
~ comparison_allTogether(reference = fadul1.2_processed,
target = fadul1.1_processed,
theta = .))
kmComparison_allCMCs_rev <- kmComparisonFeatures_rev %>%
mutate(originalMethodClassif = decision_CMC(cellIndex = cellIndex,
x = x,
y = y,
theta = theta,
corr = pairwiseCompCor,
xThresh = 20,
thetaThresh = 6,
corrThresh = .5),
highCMCClassif = decision_CMC(cellIndex = cellIndex,
x = x,
y = y,
theta = theta,
corr = pairwiseCompCor,
xThresh = 20,
thetaThresh = 6,
corrThresh = .5,
tau = 1))
```
The logic required to combine the results in `kmComparison_allCMCs` and `kmComparison_allCMCs_rev` can get a little complicated (although entirely doable using `dplyr` and other `tidyverse` functions). The user must decide precisely how results from both directions are to be combined. For example, if one direction fails the High CMC criterion yet the other passes, should we treat this as if *both* directions failed? Will you only count the CMCs in the direction that passed? We have found the best option to be treating a failure in one direction as a failure in both directions -- such cartridge case pairs would then be assigned the minimum of the two CMC counts determined under the original method of [Song (2013)](https://tsapps.nist.gov/publication/get_pdf.cfm?pub_id=911193). The `decision_combineDirections` function implements the logic to combine the two sets of results, assuming the data frame contains columns named `originalMethodClassif` and `highCMCClassif` (as defined above).
```{r}
decision_combineDirections(kmComparison_allCMCs %>%
select(-c(cellHeightValues,alignedTargetCell)),
kmComparison_allCMCs_rev)
```
The final step is to decide whether the number of CMCs computed under your preferred method is large enough to declare the cartridge case a match. [Song (2013)](https://tsapps.nist.gov/publication/get_pdf.cfm?pub_id=911193) originally proposed using a CMC count equal to 6 as the decision boundary (i.e., classify "match" if CMC count is greater than or equal to 6). This has been shown to not generalize to other proposed methods and data sets (see, e.g., [Chen et al. (2017)](https://www.sciencedirect.com/science/article/pii/S0379073817303420?via%3Dihub)). A more principled approach to choosing the CMC count has not yet been described.
Finally, we can visualize the regions of the scan identified as CMCs.
```{r}
cmcPlot(reference = fadul1.1_processed,
target = fadul1.2_processed,
cmcClassifs = kmComparisonFeatures %>%
mutate(originalMethodClassif = decision_CMC(cellIndex = cellIndex,
x = x,
y = y,
theta = theta,
corr = pairwiseCompCor,
xThresh = 20,
thetaThresh = 6,
corrThresh = .5)) %>%
group_by(cellIndex) %>%
filter(pairwiseCompCor == max(pairwiseCompCor)),
cmcCol = "originalMethodClassif")
```