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
|
# -*- coding: utf-8 -*-
#
# Identify the abnormal pixels (with relatively larger values) in the
# given X-ray count image, and replace their values with random Poisson
# values whose parameter lambda is determined by the neighboring pixels.
#
#
# Aaron LI
# 2015/09/01
# Updated: 2015/09/02
#
# Fill a vector of a row/column of X-ray image with Poisson distribution.
#
# The parameter lambda of the Poisson distribution is determined by the
# mean value of the left n (default: 3) elements.
#
# If the value of one element is greater than or equal to (>=)
# qpois(prob, lambda), then its value is replaced with rpois(1, lambda).
#
# Arguments:
# vec - input data vector
# n - number of elements used to calculate the lambda (default: 3)
# prob - quantile probability (default: 95%)
#
# Return:
# a vector of the same length with abnormal values replaced
fill.pois.vec <- function(vec, n=3, prob=0.95) {
# Select the given index of element from the vector
# Arguments:
# bc - boundary condition:
# + "cyclic": vec[0, -1] := vec[length(vec), length(vec)-1]
# + "symmetric": vec[0, -1] := vec[1, 2]
elem <- function(vec, index, bc="cyclic") {
if (index <= 0) {
if (bc == "cyclic") {
index <- length(vec) + index
} else if (bc == "symmetric") {
index <- 1 - index
} else {
stop(paste("Invalid boundary condition:", bc, "\n"))
}
}
return(vec[index])
}
# Calculate the mean value of left n elements for the given element.
mean.left <- function(vec, index, n=3, bc="cyclic") {
elements <- get(class(vec))(n)
for (i in 1:n) {
elements[i] <- elem(vec, index-i)
}
return(mean(elements))
}
#
vec.filled <- vec
for (i in 1:length(vec)) {
lambda <- mean.left(vec, i, n)
if (vec[i] >= qpois(prob, lambda)) {
vec.filled[i] <- rpois(1, lambda)
}
}
return(vec.filled)
}
# Fill a matrix of the X-ray image with Poisson distribution by row or column.
#
# For more details, see 'fill.pois.vec()'.
#
# Arguments:
# mat - input image matrix
# byrow - where Poisson fill the matrix by row (default by column)
# n - number of elements used to calculate the lambda (default: 3)
# prob - quantile probability (default: 95%)
#
# Return:
# a matrix of the same size with abnormal values replaced
fill.pois.mat <- function(mat, byrow=FALSE, n=3, prob=0.95) {
mat.filled <- mat
if (byrow) {
# process by row
rows <- nrow(mat)
for (r in 1:rows) {
vec <- mat[r, ]
vec.filled <- fill.pois.vec(vec, n, prob)
mat.filled[r, ] <- vec.filled
}
} else {
# process by column
cols <- ncol(mat)
for (c in 1:cols) {
vec <- mat[, c]
vec.filled <- fill.pois.vec(vec, n, prob)
mat.filled[, c] <- vec.filled
}
}
return(mat.filled)
}
# Identify the abnormal pixels (with relatively larger values) in the
# given X-ray count image, and replace their values with random Poisson
# values whose parameter lambda is determined by the neighboring pixels.
#
# The refilled image is the average of the two images, which are the
# original image processed with 'fill.pois.mat()' by row and column,
# respectively.
#
# TODO: to verify???
# The two-step procedure is employed to avoid the grid-like pattern/structure
# in the refilled image.
#
# For more details, see 'fill.pois.vec()'.
#
# Arguments:
# img - input image (a matrix)
# n - number of elements used to calculate the lambda (default: 3)
# prob - quantile probability (default: 95%)
#
# Return:
# a matrix of the same size with abnormal values replaced
fill.pois.img <- function(img, n=3, prob=0.95) {
img.fillbycol <- fill.pois.mat(img, byrow=FALSE, n=n, prob=prob)
img.fillbyrow <- fill.pois.mat(img, byrow=TRUE, n=n, prob=prob)
img.filled <- (img.fillbycol + img.fillbyrow) / 2
return(img.filled)
}
# vim: set ts=8 sw=4 tw=0 fenc=utf-8 ft=r: #
|
