/
core-out-sfourier.R
executable file
·232 lines (224 loc) · 8.16 KB
/
core-out-sfourier.R
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
##### Core function for radii variation Fourier analyses
#' Radii variation Fourier transform (equally spaced curvilinear abscissa)
#'
#' \code{sfourier} computes radii variation Fourier analysis from a matrix or a
#' list of coordinates where points are equally spaced aong the curvilinear abscissa.
#'
#' @param x A \code{list} or \code{matrix} of coordinates or an \code{Out} object
#' @param nb.h \code{integer}. The number of harmonics to use. If missing, 12 is used on shapes;
#' 99 percent of harmonic power on Out objects, both with messages.
#' @return A list with following components:
#' \itemize{
#' \item \code{an} vector of \eqn{a_{1->n}} harmonic coefficients
#' \item \code{bn} vector of \eqn{b_{1->n}} harmonic coefficients
#' \item \code{ao} ao harmonic coefficient
#' \item \code{r} vector of radii lengths
#' }
#' @family sfourier
#' @note The implementation is still quite experimental (as of Dec. 2016)
#' @references Renaud S, Michaux JR (2003): Adaptive latitudinal trends in the mandible shape
#' of \emph{Apodemus} wood mice. \emph{J Biogeogr} 30:1617-1628.
#' @examples
#' molars[4] %>%
#' coo_center %>% coo_scale %>% coo_interpolate(1080) %>%
#' coo_slidedirection("right") %>%
#' coo_sample(360) %T>% coo_plot(zoom=2) %>%
#' sfourier(16) %>%
#' sfourier_i() %>%
#' coo_draw(bor="red", points=TRUE)
#' @rdname sfourier
#' @export
sfourier <- function(x, nb.h){
UseMethod("sfourier")
}
#' @rdname sfourier
#' @export
sfourier.default <- function(x, nb.h){
shp <- x
# deduces the number of points
N <- coo_nb(shp)
# if missing, nb.h will be set to max (calibrate this)
if (missing(nb.h))
nb.h <- floor(N/2)
# calculates radii lengths
r <- coo_centdist(shp)
# extract an and bn coefficients using a discrete fft
an <- (Re(fft(r))/N)[1:(nb.h+1)]
bn <- (Im(fft(r))/N)[1:(nb.h+1)]
# grabs size and normalize with it
ao <- an[1]
# an[-1] <- an[-1]/ao
an <- an[-1]/ao
bn <- bn[-1]/ao
# return results
res <- list(an=an,
bn=bn,
ao=ao,
r=r)
# removes remaining names
res <- lapply(res, function(x) {names(x) <- NULL; x})
return(res)
}
#' @rdname sfourier
#' @export
sfourier.Out <- function(x, nb.h){
Out <- x
# verify
Out %<>% verify()
q <- floor(min(sapply(Out$coo, nrow)/2))
if (missing(nb.h)) {
# nb.h <- ifelse(q >= 32, 32, q)
nb.h <- calibrate_harmonicpower_sfourier(Out,
thresh = 99, plot=FALSE)$minh
# if (.is_verbose()) message("'nb.h' not provided and set to ", nb.h, " (99% harmonic power)")
}
if (nb.h > q) {
nb.h <- q # should not be 1 #todo
message("at least one outline has no more than ", q * 2,
" coordinates. 'nb.h' has been set to ", q,
" harmonics")
}
coo <- Out$coo
col.n <- paste0(rep(LETTERS[1:2], each = nb.h), rep(1:nb.h,
times = 2))
coe <- matrix(ncol = 2 * nb.h, nrow = length(coo), dimnames = list(names(coo),
col.n))
ao <- vector("numeric", length=length(coo))
for (i in seq(along = coo)) {
rf <- sfourier(coo[[i]], nb.h = nb.h) #todo: vectorize
coe[i, ] <- c(rf$an, rf$bn)
ao[i] <- rf$ao
}
res <- OutCoe(coe = coe, fac = Out$fac, method = "sfourier", norm = TRUE)
res$cuts <- ncol(res$coe)
res$ao <- ao
return(res)
}
#' @rdname sfourier
#' @export
sfourier.list <- function(x, nb.h){
lapply(x, sfourier, nb.h)
}
#' Inverse radii variation Fourier transform
#'
#' \code{sfourier_i} uses the inverse radii variation (equally spaced curvilinear abscissa) transformation to
#' calculate a shape, when given a list with Fourier coefficients, typically
#' obtained computed with \link{sfourier}.
#'
#'
#' @param rf A \code{list} with \code{ao}, \code{an} and \code{bn} components,
#' typically as returned by \code{sfourier}.
#' @param nb.h \code{integer}. The number of harmonics to calculate/use.
#' @param nb.pts \code{integer}. The number of points to calculate.
#' @param dtheta \code{logical}. Whether to use the dtheta correction method.
#' \code{FALSE} by default. When \code{TRUE}, tries to correct the angular difference between
#' reconstructed points; otherwise equal angles are used.
#' @return A list with components: \item{x }{\code{vector} of
#' \code{x}-coordinates.} \item{y }{\code{vector} of \code{y}-coordinates.}
#' \item{angle}{\code{vector} of angles used.} \item{r}{\code{vector} of radii
#' calculated.}
#' @family sfourier
#' @references Renaud S, Pale JRM, Michaux JR (2003): Adaptive latitudinal trends in the mandible shape
#' of \emph{Apodemus} wood mice. \emph{J Biogeogr} 30:1617-1628.
#' @examples
#' coo <- coo_center(bot[1]) # centering is almost mandatory for sfourier family
#' coo_plot(coo)
#' rf <- sfourier(coo, 12)
#' rf
#' rfi <- sfourier_i(rf)
#' coo_draw(rfi, border='red', col=NA)
#'
#' @export
sfourier_i <-
function(rf, nb.h, nb.pts = 120, dtheta=FALSE) {
if (!all(c("an", "bn") %in% names(rf))) {
stop("a list containing 'an' and 'bn' harmonic coefficients must be provided")
}
ao <- ifelse(is.null(rf$ao), 1, rf$ao)
an <- rf$an
bn <- rf$bn
if (missing(nb.h)) {
nb.h <- length(an)
}
if (nb.h > length(an)) {
nb.h <- length(an)
message("nb.h cannot be higher than length(rf$an) and has been set to ", nb.h)
}
theta <- seq(0, 2 * pi, length = nb.pts)
harm <- matrix(NA, nrow = nb.h, ncol = nb.pts)
for (i in 1:nb.h) {
harm[i, ] <- an[i] * cos(i * theta) + bn[i] * sin(i * theta)
}
r <- ao/2 + apply(harm, 2, sum)
# theta corrected using ds/dtheta method
if (dtheta){
ds <- (2*pi)/sum(1/r)
theta <- numeric(length(r))
for (i in 2:nb.pts){
theta[i] <- theta[i-1] + ds/r[i]
}
coo <- cbind(r*cos(theta), r*sin(theta))
} else {
# theta based on the circle
Z <- complex(modulus = r, argument = theta)
coo <- cbind(Re(Z), Im(Z))
}
colnames(coo) <- c("x", "y")
return(coo)
}
#' Calculates and draw 'sfourier' shapes.
#'
#' \code{sfourier_shape} calculates a 'Fourier radii variation shape' given
#' Fourier coefficients (see \code{Details}) or can generate some 'sfourier'
#' shapes.
#'
#' \code{sfourier_shape} can be used by specifying \code{nb.h} and
#' \code{alpha}. The coefficients are then sampled in an uniform distribution
#' \eqn{(-\pi ; \pi)} and this amplitude is then divided by
#' \eqn{harmonicrank^alpha}. If \code{alpha} is lower than 1, consecutive
#' coefficients will thus increase. See \link{sfourier} for the mathematical
#' background.
#'
#' @param an \code{numeric}. The \eqn{a_n} Fourier coefficients on which to
#' calculate a shape.
#' @param bn \code{numeric}. The \eqn{b_n} Fourier coefficients on which to
#' calculate a shape.
#' @param nb.h \code{integer}. The number of harmonics to use.
#' @param nb.pts \code{integer}. The number of points to calculate.
#' @param alpha \code{numeric}. The power coefficient associated with the
#' (usually decreasing) amplitude of the Fourier coefficients (see
#' \bold{Details}).
#' @param plot \code{logical}. Whether to plot or not the shape.
#' @return A matrix of (x; y) coordinates.
#' @family sfourier
#' @references Renaud S, Pale JRM, Michaux JR (2003): Adaptive latitudinal trends in the mandible shape
#' of \emph{Apodemus} wood mice. \emph{J Biogeogr} 30:1617-1628.
#' @examples
#' rf <- sfourier(bot[1], 24)
#' sfourier_shape(rf$an, rf$bn) # equivalent to sfourier_i(rf)
#' sfourier_shape() # not very interesting
#'
#' sfourier_shape(nb.h=12) # better
#' sfourier_shape(nb.h=6, alpha=0.4, nb.pts=500)
#'
#' # Butterflies of the vignette' cover
#' panel(Out(a2l(replicate(100,
#' sfourier_shape(nb.h=6, alpha=0.4, nb.pts=200, plot=FALSE)))))
#' @export
sfourier_shape <- function(an, bn, nb.h, nb.pts = 80, alpha = 2,
plot = TRUE) {
if (missing(nb.h) & missing(an))
nb.h <- 6
if (missing(nb.h) & !missing(an))
nb.h <- length(an)
if (missing(an))
an <- runif(nb.h, -pi, pi)/(1:nb.h)^alpha
if (missing(bn))
bn <- runif(nb.h, -pi, pi)/(1:nb.h)^alpha
rf <- list(an = an, bn = bn, ao = 0)
shp <- sfourier_i(rf, nb.h = nb.h, nb.pts = nb.pts)
if (plot)
coo_plot(shp)
return(shp)
}
##### end sfourier