I’m considering implementing Kikuchi line drawing for my studies at https://github.com/SaitouToshihide/rgl.cry/tree/DiffractionCone.
dp_demo(system.file("cubic_i.cif", package = "rgl.cry"), zoom = 0.03, res = 1.3)
##rgl::snapshot3d("cubic_i_01.png", fmt="png", width=500, height=500, webshot = FALSE)
dp_demo(system.file("cubic_i.cif", package = "rgl.cry"), zoom = 0.0005, res = 1.3)
##rgl::snapshot3d("cubic_i_02.png", fmt="png", width=500, height=500, webshot = FALSE)
library(devtools)
install_github("SaitouToshihide/rgl.cry", ref="5e63cad")
Run dp_demo()
.
dp_demo() # default CIF data are loaded.
Let’s see another crystal system.
dp_demo(system.file("triclinic_p.cif", package = "rgl.cry"), zoom = 0.03, res = 2.0)
##rgl::snapshot3d("triclinic_p_001_01.png", fmt="png", width=500, height=500, webshot = FALSE)
dp_demo(system.file("triclinic_p.cif", package = "rgl.cry"), zoom = 0.0003, res = 2.0)
##rgl::snapshot3d("triclinic_p_001_02.png", fmt="png", width=500, height=500, webshot = FALSE)
In this case, even in a triclinic lattice, diffraction spots such as (1 1 0) are visible. This is because the corresponding [1 1 0] vector in reciprocal space remains parallel to the xy-plane, even if the real-space [1 1 0] vector is inclined relative to the xy-plane. We should not focus on imaging the vector denoted by Miller indices in real space. Instead, we should consider the repetition of planes and their normal direction or alternatively, image a vector in reciprocal space denoted by Miller indices. Let’s see.
cry_demo()
align("90 0")
align("0 50") #
library(pryr)
lCIF <- rgl.cry:::pkg$inst[[nrow(rgl.cry:::pkg$inst), "lCIF"]]
a <- lCIF$HEADER$CELL$A$VAL
b <- lCIF$HEADER$CELL$B$VAL
c <- lCIF$HEADER$CELL$C$VAL
aa <- lCIF$HEADER$CELL$ALPHA$VAL # degree
bb <- lCIF$HEADER$CELL$BETA$VAL
cc <- lCIF$HEADER$CELL$GAMMA$VAL
mat02 <- cry::xtal_mat02(a, b, c, aa, bb, cc)
vertices <- rbind(c(1,0,0), c(1,0,1), c(0,1,1), c(0,1,0))
vertices <- t(apply(vertices, 1, function(v) v %*% mat02))
combinations <- c(1,2,3,1, 1,4,3,1)
rgl::shade3d(rgl::mesh3d(vertices, quads = combinations), col = "blue", alpha = 0.9)
vertices <- rbind(c(1,0,0), c(1,0,1), c(0,1,1), c(0,1,0))
vertices <- t(apply(vertices, 1, function(v) v + c(1, 1, 0)))
vertices <- t(apply(vertices, 1, function(v) v %*% mat02))
rgl::shade3d(rgl::mesh3d(vertices, quads = combinations), col = "green", alpha = 0.9)
##rgl::snapshot3d("triclinic_p.png", fmt="png", width=500, height=500, webshot = FALSE)
This is the side view of the image looking down along with c-axis indicated above as dp pattern. We can observe that the direction [1 1 0] inclined in real space. However, as the repetition of planes, this lies parallel to the xy-plane.
IMPORTANT: Adding 3D graphics objects is made easy and flexible within R using the rgl and cry packages.
Let’s check this numerically.
> mat02 <- cry::xtal_mat02(a, b, c, aa, bb, cc)
> ea1 <- as.numeric(c(1, 0, 0) %*% mat02)
> ea2 <- as.numeric(c(0, 1, 0) %*% mat02)
> ea3 <- as.numeric(c(0, 0, 1) %*% mat02)
> library(pracma)
> V <- as.numeric(ea1 %*% cross(ea2, ea3))
> eb1 <- cross(ea2, ea3) / V
> eb2 <- cross(ea3, ea1) / V
> eb3 <- cross(ea1, ea2) / V
## The [1 1 0] vector in real space and its corresponding coordinates in the Cartesian system.
> 1*ea1 + 1*ea2 + 0*ea3
[1] 6.554986 9.455095 5.356230
## The [1 1 0] vector in reciprocal space and its corresponding coordinates in the Cartesian system.
> 1*eb1 + 1*eb2 + 0*eb3
[1] 0.1020649 0.1407670 0.0000000
In reciprocal space, the vector denoted by Miller indices is perpendicular to the plane in real space with the same Miller indices.
You may be confused if you are asked to calculate the interplanar spacing in a triclinic lattice. Don’t worry, it’s easy when using reciprocal lattice vectors.
## The interplanar spacing of the (1 1 0) planes in this crystal is 1/0.1738753 Å.
> hkl[which(hkl$H=="1" & hkl$K=="1" & hkl$L=="0"),]
H K L S
124 1 1 0 0.1738753
> 1 / 0.1738753
[1] 5.751248
> n110b <- (1*eb1 + 1*eb2 + 0*eb3)
> norm(n110b, "2")
[1] 0.1738753
## On the other hand,
> n110a <- (1*ea1 + 1*ea2 + 0*ea3) * c(1, 1, 0)
> norm(n110a, "2")
[1] 11.50507
> norm(n110a, "2") / 2
[1] 5.752536
## Although the calculated value may initially appear plausible after being halved,
## this approach, using the lattice vector of real space, is fundamentally flawed and unproductive.
Prepare image.
library(rgl.cry)
dp_demo(system.file("cubic_i.cif", package = "rgl.cry"), zoom = 0.03, res = 1.3)
## Remove spot, label, etc.
rgl::pop3d(tag = c("rlpoints0", "rlpoints1", "rlpoints2"))
objIds <- rgl::ids3d(tag = TRUE, subscene = 0)
objIds <- objIds[grep("^rlpoints[4-5].*", objIds$tag), ]
rgl::pop3d(id = objIds$id)
## remove manually in this way
rgl::ids3d(subscene = 0)
##rgl::pop3d(id = 16)
rgl::par3d(zoom = 0.4)
## Save the image.
rgl::snapshot3d("cubic_i_03.png", fmt="png", width=500, height=500, webshot = FALSE)
Utilize ImageMagick to blur, grayscale conversion, and inversion to achieve a visual representation that resembles the original observed image. ((In this case, procedures based on theory can be considered, but in general image applications, what methods are there for converting basic features into actual images, and what would happen if they were applied? Applying image blur can be motivated by our visual perception or by limitations inherent in the observation process, making it a common image effect.))
% convert -blur 10x10 cubic_i_03.png -colorspace Gray -negate cubic_i_03bw.png
You can also use the magick package.
library(magick)
img <- image_read("cubic_i_03.png")
img <- image_blur(img, 10, 10)
img <- image_convert(img, type = 'grayscale')
img <- image_negate(img)
image_browse(img)
image_write(img, path = "cubic_i_03bw.png", format = "png")
Reconstruct the line from Hough transform data and overlay it onto the original image.
library(imager)
file <- "cubic_i_03bw.png"
im <- load.image(file)
plot(im)
hough_line(im, ntheta=500) %>% plot # Very hard to see
df <- hough_line(im, ntheta=500, data.frame=TRUE)
df <- subset(df, score > quantile(score, .9991))
df <- df[which(df$score != 0),]
##png(filename="cubic_i_03ovl.png", width=500, height=500, units="px")
plot(im)
nfline(df$theta, df$rho, col="blue")
##dev.off()
To improve the visual representation of Hough transform data, I used the following settings:
library(ggplot2)
##png(filename="cubic_i_03ht.png", width=500, height=500, units="px")
ggplot(data = df) +
geom_point(aes(x = theta, y = rho, color = score), size = 0.1) +
scale_color_gradient(low = "#56B1F7", high = "#132B43") +
theme(panel.background = element_rect(fill = "gray90"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.background = element_rect(fill = "white"))
##dev.off()
We detect lines in the image using Hough transform and keep only the highest-scoring 1-0.9991 %. Let’s check this.
> plot(im)
> nfline(df$theta, df$rho, col="blue")
## To plot the detected line, substitute the values of theta and rho obtained
## from the Hough transform into the equation cos(theta)*x + sin(theta)*y = rho.
## This can be done using nfline. Let's check it out.
> df[which(df$score == max(df$score)), ]
theta rho score
442140 3.928565 -362.169 11.64875
>
## y = (rho - cos(theta)*x)/sin(theta)
## y = (rho - cos(3.928565)*x)/sin(3.928565)
> f1 <- function(x){ (-362.169 + 0.7059928*x) / -0.708219 }
> par(new=T)
> plot(f1, n=1000, col="red", xaxt="n", yaxt="n", xlab="", ylab="")
Utilize the pracma package for peak detection.
library(pracma)
im <- load.image(file)
df <- hough_line(im, ntheta=500, data.frame=TRUE)
## Treating a vertical raster column (rho) as a one-dimensional data.
## Then find peaks using findpeaks (adjust minpeakheight and threshold as needed).
data_vector_rho <- df$score
peak_indices_rho <- findpeaks(data_vector_rho, minpeakheight = 2.5, threshold = 0.5)
df4 <- df[peak_indices_rho[,2],]
##png(filename="cubic_i_04ovl.png", width=500, height=500, units="px")
plot(im, xaxs = "i", yaxs = "i")
nfline(df4$theta, df4$rho, col="red")
##dev.off()
##png(filename="cubic_i_04ht.png", width=500, height=500, units="px")
ggplot(data = df4) +
geom_point(aes(x = theta, y = rho, color = score), size = 0.1) +
scale_color_gradient(low = "#56B1F7", high = "#132B43") +
theme(panel.background = element_rect(fill = "gray90"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.background = element_rect(fill = "white"))
##dev.off()
The four major lines (\, /, |, -) in the image on the left share identical slope and Y-intercept. This characteristic corresponds to the existence of 4 major spots in the image on the right.
And this does not lead to further improvement:
## Treating a horizontal raster row (theta).
data_vector_theta <- df[order(df[, "rho"], df[, "theta"]), "score"]
peak_indices_theta <- findpeaks(data_vector_rho, minpeakheight = 2.5, threshold = 0.5)
nrho <- nrow(df[which(df$theta == 0),])
row_coords <- ceiling(peak_indices_theta[,2] / nrho)
col_coords <- peak_indices_theta[,2] %% nrho
df5 <- df[nrho * row_coords + col_coords,]
## Or even further, multiply by the vertical result?