在 Photoshop 或 R 中将来自耦合显微镜/光谱学的图像和数据拼接成全景图
Cai*_*lin 18 javascript photoshop r raster extendscript
我有一组由耦合扫描电子显微镜和能量色散光谱生成的图像和 X 射线数据。这是我的问题：
我像这样对岩石表面的横断面进行了成像（紫色框勾勒出横断面区域）：
我想要真正的高分辨率，所以我使用 7 张放大 3000 倍的图像并将它们与 Photoshop 中的 photomerge 脚本拼接在一起。这是单个图像的示例：
以及它在光合图像横断面上的位置：
在这 7 个位置中的每一个，我还收集了 X 射线数据，为检测到的每个元素生成元素图，并将其写入 TIFF。我还想将每个元素图 TIFF 拼接在一起，以便我可以将其叠加在岩石的合并横断面图像上。这是我想要的结果：
问题是元素贴图没有足够的功能，无法使用 photomerge 将它们拼接在一起。它基本上是二进制的 - 如果检测到元素，则像素是某种颜色（例如，在我的示例图像中，红色代表铁，黄色代表硫），如果未检测到元素，则像素为黑色。您可以看到大部分元素地图大部分是黑色的。
我现在有 ~20 个横断面 x 7 个图像，每个 x ~10 个元素，这导致需要将 ~1400 个图像放在一起，因此需要自动化。
我的想法是用 photomerge 将岩石图像拼接在一起。photomerge 的输出是一个智能对象，其中每个图像都是一个图层。然后，我将使用一个脚本来获取 photomerged 图像对象中 7 个图像中每一个的左上角坐标、宽度和高度。然后，我会将这些属性放置并分配给 7 个图像的每个相应元素映射，以生成“合并”元素映射以覆盖在图像上。我试图自己解决这个问题，但我不精通 javascript 并且无法围绕 Photoshop API。
我在这里在 Github 上上传了一个示例数据集。7 个横断面位置从左到右：-2、-1、0、1、2、3、4。每个位置都有岩石和子目录的图像以及元素数据。
我不知道 Photoshop 或 R，JavaScript 也是如此：
const names = { // map from directory names to patterns (where "#" stands for position index) of names of images therein
"SEM_images" : "pos# image.tif",
"Al" : "Al K?1 pos# map data.tif",
"Ba" : "Ba L?1 pos# map data.tif",
"C" : "C K?1_2 pos# map data.tif",
"Ca" : "Ca K?1 pos# map data.tif",
"Fe" : "Fe K?1 pos# map data.tif",
"Hg" : "Hg L?1 pos# map data.tif",
"Ir" : "Ir L?1 pos# map data.tif",
"K" : "K K?1 pos# map data.tif",
"Mg" : "Mg K?1_2 pos# map data.tif",
"Mn" : "Mn K?1 pos# map data.tif",
"Na" : "Na K?1_2 pos# map data.tif",
"O" : "O K?1 pos# map data.tif",
"Os" : "Os L?1 pos# map data.tif",
"P" : "P K?1 pos# map data.tif",
"S" : "S K?1 pos# map data.tif",
"Si" : "Si K?1 pos# map data.tif",
"Ti" : "Ti K?1 pos# map data.tif"
}
const SCALE = 1/10 // scale of output images
const OVERLAP = 1.0 // minimum *tested* (horizontal) overlap of images relative to their width
const H_BOX = 0.1 // height of comparison box relative to height of images
const W_BOX = 0.1 // width of comparison box relative to width of images
const ADJUSTMENT = 0 // (vertical) adjustment of comparison box [pixels]
/* Merge images given:
* dataset - dataset address as String
* directory - directory name for images as String
* pattern - pattern (where "#" stands for position index) of names of images in directory
* pos_min - minimum position index of images as Number
* pos_max - maximum position index of images as Number
*/
function Merge(dataset, directory, pos_min, pos_max) {
if (dataset[dataset.length - 1] != "/") dataset += "/"
const images = []
for (let pos = pos_min; pos <= pos_max; ++pos) (images[pos - pos_min] = new Image).src = dataset + directory + "/" + names[directory].replace("#", pos)
merge(images, dataset, pos_min, pos_max)
}
function Laplacian(imagedata) { // 5-point stencil approximation
const data = imagedata.data
const L = data.length/4
const grayscale = new Float32Array(L)
for (let i = 0; i < L; ++i) {
const I = 4*i
grayscale[i] = (data[I ] + data[I + 1] + data[I + 2])/3
}
const Laplacian = new Float32Array(L)
//const H = imagedata.height
const Hm1 = imagedata.height - 1
const W = imagedata.width
const Wm1 = W - 1
for (let r = 1; r < Hm1; ++r) {
const R = r*W
for (let c = 1; c < Wm1; ++c) {
const i = R + c
Laplacian[i] = grayscale[i - W] + grayscale[i + W] + grayscale[i - 1] + grayscale[i + 1] - 4*grayscale[i]
}
}
for (let c = 1; c < Wm1; ++c) {
//const i = c
Laplacian[c] = grayscale[c + W] + grayscale[c - 1] + grayscale[c + 1] - 4*grayscale[c]
}
for (let c = 1; c < Wm1; ++c) {
const i = Hm1*W + c
Laplacian[i] = grayscale[i - W] + grayscale[i - 1] + grayscale[i + 1] - 4*grayscale[i]
}
for (let r = 1; r < Hm1; ++r) {
const i = r*W
Laplacian[i] = grayscale[i - W] + grayscale[i + W] + grayscale[i + 1] - 4*grayscale[i]
}
for (let r = 1; r < Hm1; ++r) {
const i = r*W + Wm1
Laplacian[i] = grayscale[i - W] + grayscale[i + W] + grayscale[i - 1] - 4*grayscale[i]
}
{
const Lm1 = L - 1
const LmW = L - W
Laplacian[0 ] = grayscale[W ] + grayscale[1 ] - 4*grayscale[0 ]
Laplacian[W ] = grayscale[2*W ] + grayscale[Wm1 ] - 4*grayscale[W ]
Laplacian[LmW] = grayscale[LmW - W] + grayscale[LmW + 1] - 4*grayscale[LmW]
Laplacian[Lm1] = grayscale[Lm1 - W] + grayscale[Lm1 - 1] - 4*grayscale[Lm1]
}
return Laplacian
}
function merge(images, dataset, pos_min, pos_max) {
for (const image of images) if (!image.complete) {
setTimeout(merge, 1000, images, dataset, pos_min, pos_max) // wait 1000ms = 1s
return
}
let Row, Col
const Coords = [[Row = 0, Col = 0]]
let index = 0
let image = images[index]
const H = image.naturalHeight
const W = image.naturalWidth
if (W*H == 0) return []
const canvas = document.createElement("canvas")
canvas.height = H
canvas.width = W
const context = canvas.getContext('2d')
context.drawImage(image, 0, 0)
let prev = Laplacian(context.getImageData(0, 0, W, H))
const length = images.length
const h = Math.round(H_BOX*H)
const Hmh = H - h
const w = Math.round(W_BOX*W)
const o = Math.max(Math.round((1 - OVERLAP)*W), w)
const Wmw = W - w
const row_offset = Math.round(Hmh/2) + ADJUSTMENT
const offset = row_offset*W
for (++index; index < length; ++index) {
image = images[index]
if (image.naturalHeight != H || image.naturalWidth != W) alert("Dimension mismatch: " + image.src)
context.drawImage(image, 0, 0)
const curr = Laplacian(context.getImageData(0, 0, W, H))
let max = -1
let row, col
for (let r = 0; r < Hmh; ++r) {
const R = r*W
for (let c = o; c < Wmw; ++c) {
let m = 0
for (let i = 0; i < h; ++i) {
const I = i*W
const K = R + I + c
const k = offset + I
for (let j = 0; j < w; ++j) if (prev[K + j]*curr[k + j] > 0) ++m
}
if (m > max) {
max = m
row = r
col = c
}
}
}
Coords[index] = [(Row += row - row_offset)/H, (Col += col)/W]
prev = curr
}
Stitch(dataset, pos_min, pos_max, Coords)
}
function Stitch(dataset, pos_min, pos_max, Coords) {
if (dataset[dataset.length - 1] != "/") dataset += "/"
document.body.appendChild(document.createElement("h1")).innerText = `${dataset} :[${pos_min},${pos_max}] @${JSON.stringify(Coords)}`
const tasks = []
for (const directory in names) {
document.body.appendChild(document.createElement("h2")).innerText = directory
const images = []
for (let pos = pos_min; pos <= pos_max; ++pos) (images[pos - pos_min] = new Image).src = dataset + directory + "/" + names[directory].replace("#", pos)
const target = document.body.appendChild(document.createElement("img"))
target.height = 0
target.width = 0
tasks.push([images, target])
}
process(tasks, Coords)
}
const ROW = 0
const COL = 1
function stitch(images, Coords) {
let image
let index
for (index in images) {
image = images[index]
if (image.naturalHeight != 0 && image.naturalWidth != 0) break
}
const H = image.naturalHeight
const W = image.naturalWidth
const canvas = document.createElement("canvas")
let r_min = 0
let r_max = 0
let c_min = 0
let c_max = 0
for (coords of Coords) {
const r = coords[ROW]
const c = coords[COL]
if (r < r_min) r_min = r
if (r > r_max) r_max = r
if (c < c_min) c_min = c
if (c > c_max) c_max = c
}
canvas.height = (r_max - r_min + 1)*H
canvas.width = (c_max - c_min + 1)*W
const context = canvas.getContext('2d')
if (context == null) {
let list = ""
for (const image of images) list += "\n- " + image.src
alert("Too large: stitching area required for:" + list)
return
}
let coords = Coords[index]
let row = (coords[ROW] - r_min)*H
let col = (coords[COL] - c_min)*W
context.drawImage(image, col, row)
const length = images.length
for (++index; index < length; ++index) {
image = images[index]
if (image.naturalHeight == 0 || image.naturalWidth == 0) continue
if (image.naturalHeight != H || image.naturalWidth != W) alert("Dimension mismatch: " + image.src)
coords = Coords[index]
row = coords[ROW]*H
col = coords[COL]*W
context.drawImage(image, col, row)
}
return canvas.toDataURL()
}
function process(tasks, Coords) {
const task = tasks.shift()
const images = task[0]
for (const image of images) if (!image.complete) {
tasks.push(task)
setTimeout(process, 1000, tasks, Coords) // wait 1000ms = 1s
return
}
const target = task[1]
target.src = stitch(images, Coords)
target.onload = function () {
this.height = SCALE*this.naturalHeight
this.width = SCALE*this.naturalWidth
this.style = "border: solid black 1px"
}
if (tasks.length > 0) process(tasks, Coords)
}
要运行，请执行以下操作：
Merge("https://raw.githubusercontent.com/CaitlinCasar/dataStitcher/master/example_dataset/", "SEM_images", -2, 4)
带有 Fe 覆盖层的 SEM_images 示例：
如果完全在 Photoshop 脚本中或完全在 R 中完成此操作，那就太好了，但我有一个同时使用这两种脚本的可行解决方案。我手动从 Photoshop 的 photomerge 输出中的每个图层的左下角拉出 XY 坐标，并将它们分配给 R 中的光栅化图像：
\n#!/usr/bin/env Rscript\n\n#Caitlin Casar\n#github.com/caitlincasar\n#19 March 2020\n#DataStitchR stitches panoramic images of SEM images coupled to x-ray energy dispersive spectroscopy.\n\n# usage: ./dataStitchR. -f "/Users/Caitlin/Desktop/dataStitcher/example_dataset" -b "/Users/Caitlin/Desktop/dataStitcher/example_dataset/SEM_images" -c "/Users/Caitlin/Desktop/dataStitcher/coordinates.txt" -z ".tif" -m ".tif" -a "Unknown|Os|SEM" -d "overview" -y "overview" -p TRUE -o "example" -n "example"\n# usage Rscript dataStitchR.R -f "/Users/Caitlin/Desktop/dataStitcher/example_dataset" -b "/Users/Caitlin/Desktop/dataStitcher/example_dataset/SEM_images" -c "/Users/Caitlin/Desktop/dataStitcher/coordinates.txt" -z ".tif" -m ".tif" -a "Unknown|Os|SEM" -d "overview" -y "overview" -p TRUE -o "example" -n "example"\n\nmessage("\n ##### ###### \n ##### ## ##### ## # # ##### # ##### #### # # # # \n # # # # # # # # # # # # # # # # # \n # # # # # # # ##### # # # # ###### ###### \n # # ###### # ###### # # # # # # # # # \n # # # # # # # # # # # # # # # # # # \n ##### # # # # # ##### # # # #### # # # # \n \n\n ")\n\nmessage("DataStitchR was created by Caitlin Casar and is maintained at github.com/CaitlinCasar/dataStitchR.\n ")\nmessage("DataStitchR stitches panoramic images of SEM images coupled to x-ray energy dispersive spectroscopy.\n ")\nmessage("For help, run 'dataStitchR.R --help'.\n \n \n ")\n\nsuppressPackageStartupMessages(require(optparse))\n\noption_list = list(\n make_option(c("-f", "--file"), action="store", default=getwd(), type='character',\n help="Input parent directory with subdirectories of element xray data to be stitched. The element names should be abbreviated, i.e. 'Ca' for calcium."),\n make_option(c("-o", "--out"), action="store", default=NA, type='character',\n help="Output file directory. This is where your x-ray raster brick and output figures will be saved."),\n make_option(c("-n", "--name"), action="store", default=NA, type='character',\n help="Optional name for output files."),\n make_option(c("-b", "--base-images"), action="store", default=NA, type='character',\n help="SEM image file directory."),\n make_option(c("-c", "--coords"), action="store", default=NA, type='character',\n help="Tab-delimited file of xy coordinates for each image. A third column should denote stitching positions that correspond to the file names for each image."),\n make_option(c("-u", "--use-positions"), action="store", default="-?(?<![K\xce\xb11||K\xce\xb11_2])\\\\d+", type='character',\n help="Optional regex pattern to extract position IDs from each file name that corresponds to positions in the xy file. The default searches for numbers that appear after 'K\xce\xb11' or 'K\xce\xb12'. Numbers can include signs, i.e. -1 is acceptable."),\n make_option(c("-z", "--z-format"), action="store", default="*", type='character',\n help="Optional regex pattern of x-ray image formats to select for stitching, i.e. '.tif'."),\n make_option(c("-m", "--make"), action="store", default="*", type='character',\n help="Optional regex pattern of SEM image formats to select for stitching, i.e. '.tif'. You do not need to specify this unless you are generating a PDF output."),\n make_option(c("-a", "--all-exclude"), action="store", default=NA, type='character',\n help="Optional regex pattern of x-ray file directories to exclude from stitiching, i.e. the element your sample was coated with."),\n make_option(c("-d", "--drop"), action="store", default=NA, type='character',\n help="Optional regex pattern of files to exclude from x-ray data stitching."),\n make_option(c("-y", "--y-exclude"), action="store", default=NA, type='character',\n help="Optional regex pattern of files to exclude from SEM image stitiching. You do not need to specify this unless you are generating a PDF output."),\n make_option(c("-v", "--verbose"), action="store_true", default=TRUE,\n help="Print updates to console [default %default]."),\n make_option(c("-q", "--quiet"), action="store_false", dest="verbose",\n help="Do not print anything to the console."),\n make_option(c("-p", "--pdf"), action="store", default=FALSE,\n help="Generate PDF of x-ray brick colored by element superimposed on the SEM image, default is TRUE [default %default].") \n)\nopt = parse_args(OptionParser(option_list=option_list))\n\n\n#load dependencies \nsuppressPackageStartupMessages(require(pacman))\npacman::p_load(raster, magick, tidyverse, rasterVis, ggnewscale, Hmisc, cowplot)\n\n#create list all sub-directories within main directory\ndirectories <- list.dirs(opt$f, full.names = T , recursive =F)\nif(!is.na(opt$a)){\n directories <- directories[!str_detect(directories, opt$a)]\n}\n\n\n#set image coordinates\nxy <- read_delim(opt$c, delim = "\\t", col_types = cols())\npositions <- xy %>%\n select(-x, -y)\n\n\n# stitch xray images ------------------------------------------------------\n\n#stitch xray images into panoramas and store in raster brick \nxray_brick_list <- list()\nxray_data <- list()\nfor(j in 1:length(directories)){\n path = directories[j]\n files <- list.files(path, full.names = T, pattern = opt$z)\n if(!is.na(opt$d)){\n files <- files[!str_detect(files, opt$d)]\n }\n if(length(files) >0){\n xray_data[[j]] <- str_extract(path, "([^/]+$)")\n message(paste0("Stitching ", str_extract(xray_data[[j]], "([^/]+$)"), " data (element ", j, " of ", length(directories), ")..."))\n xy_id <- which(positions[[1]] %in% str_extract(files, opt$u))\n panorama <- list()\n for(i in 1:length(files)){\n message(paste0("Processing image ", i, " of ", length(files), "..."))\n image <- files[i] %>% image_read() %>% \n image_quantize(colorspace = 'gray') %>% \n image_equalize() \n temp_file <- tempfile()\n image_write(image, path = temp_file, format = 'tiff')\n image <- raster(temp_file) %>%\n cut(breaks = c(-Inf, 150, Inf)) - 1\n image <- aggregate(image, fact=4)\n image_extent <- extent(matrix(c(xy$x[xy_id[i]], xy$x[xy_id[i]] + 1024, xy$y[xy_id[i]], xy$y[xy_id[i]]+704), nrow = 2, ncol = 2, byrow = T))\n image_raster <- setExtent(raster(nrows = 704, ncols = 1024), image_extent, keepres = F)\n values(image_raster) <- values(image)\n panorama[[xy_id[i]]] <- image_raster\n }\n empty_xy_id <- which(!positions[[1]] %in% str_extract(files, opt$u))\n if(length(empty_xy_id) > 0){\n for(k in 1:length(empty_xy_id)){\n empty_raster_extent <- extent(matrix(c(xy$x[empty_xy_id[k]], xy$x[empty_xy_id[k]] + 1024, xy$y[empty_xy_id[k]], xy$y[empty_xy_id[k]]+704), nrow = 2, ncol = 2, byrow = T))\n empty_raster <- setExtent(raster(nrows = 704, ncols = 1024), empty_raster_extent, keepres = F)\n values(empty_raster) <- 0\n panorama[[empty_xy_id[k]]] <- empty_raster\n }\n }\n panorama_merged <- do.call(merge, panorama)\n xray_brick_list[[j]] <- panorama_merged\n }\n}\n\nmessage("Stitching complete. Creating x-ray brick...")\nxray_brick <- do.call(brick, xray_brick_list)\nnames(xray_brick) <- unlist(xray_data)\nmessage("...complete.")\n\n#write the brick \nmessage("Writing brick...")\n\n\nif(!is.na(opt$o)){\n dir.create(opt$o)\n if(!is.na(opt$n)){\n out_brick <- writeRaster(xray_brick, paste0(opt$o, "/", opt$n,"_brick.grd"), overwrite=TRUE, format="raster")\n x <- writeRaster(xray_brick, paste0(opt$o, "/", opt$n,"_brick.tif"), overwrite=TRUE, format="GTiff",options=c("INTERLEAVE=BAND","COMPRESS=LZW"))\n }else{\n out_brick <- writeRaster(xray_brick, path = opt$o, "brick.grd", overwrite=TRUE, format="raster")\n x <- writeRaster(xray_brick, path = opt$o, "brick.tif", overwrite=TRUE, format="GTiff",options=c("INTERLEAVE=BAND","COMPRESS=LZW"))\n }\n}else{\n if(!is.na(opt$n)){\n out_brick <- writeRaster(xray_brick, paste0(opt$n,"_brick.grd"), overwrite=TRUE, format="raster")\n x <- writeRaster(xray_brick, paste0(opt$n,"_brick.tif"), overwrite=TRUE, format="GTiff",options=c("INTERLEAVE=BAND","COMPRESS=LZW"))\n }else{\n out_brick <- writeRaster(xray_brick, "brick.grd", overwrite=TRUE, format="raster")\n x <- writeRaster(xray_brick, "brick.tif", overwrite=TRUE, format="GTiff",options=c("INTERLEAVE=BAND","COMPRESS=LZW"))\n }\n}\n\n\nmessage("...complete.")\n\n#flush everything we don't need from memory\nremove(list = c("x", "xray_brick_list", "xray_data", "empty_raster", "empty_raster_extent",\n "i", "j", "k", "path", "positions", "xy_id", \n "panorama_merged", "panorama", "empty_xy_id", "files", "directories"))\n\n\n# create base SEM image ---------------------------------------------------\nif(opt$p){\n \nif(!is.na(opt$b)){\n SEM_images <- list.files(opt$b, full.names = T, pattern = opt$m)\n}else{\n message("Please provide a file path for your SEM images to stitch. For help, see 'dataStitchR.R --help.'")\n break\n}\n\nif(!is.na(opt$y)){\n SEM_images <- SEM_images[!str_detect(SEM_images , opt$y)]\n}\nSEM_panorama <- list()\n\nmessage("Stitching SEM images into panorama...")\nfor(i in 1:length(SEM_images)){\n image <- SEM_images[i] %>% image_read() %>%\n image_quantize(colorspace = 'gray') %>%\n image_equalize()\n temp_file <- tempfile()\n image_write(image, path = temp_file, format = 'tiff')\n image <- raster(temp_file)\n image_extent <- extent(matrix(c(xy$x[i], xy$x[i] + 1024, xy$y[i], xy$y[i]+704), nrow = 2, ncol = 2, byrow = T))\n image_raster <- setExtent(raster(nrows = 704, ncols = 1024), image_extent, keepres = F)\n values(image_raster) <- values(image)\n SEM_panorama[[i]] <- image_raster\n}\nSEM_panorama_merged <- do.call(merge, SEM_panorama)\nmessage("...complete.")\n\n#flush everything we don't need from memory\nremove(list = c("SEM_panorama", "SEM_images", "image", "image_extent", "image_raster"))\n\n#optinal plot to check if SEM image merge looks correct \n#plot(SEM_panorama_merged, col = gray.colors(10, start = 0.3, end = 0.9, gamma = 2.2, alpha = NULL))\n\n\n# plot the data -----------------------------------------------------------\n\nmessage("Generating plot...")\n# Set color palette\n\n#palette source: https://sciencenotes.org/molecule-atom-colors-cpk-colors/\n\nelement_colors <- c("#FFFFFF", "#D9FFFF", "#CC80FF", "#C2FF00", "#FFB5B5", "#909090", "#3050F8",\n "#FF0D0D", "#90E050", "#B3E3F5", "#AB5CF2", "#8AFF00", "#BFA6A6", "#F0C8A0",\n "#FF8000", "#FFFF30", "#1FF01F", "#80D1E3", "#8F40D4", "#3DFF00", "#E6E6E6",\n "#BFC2C7", "#A6A6AB", "#8A99C7", "#9C7AC7", "#E06633", "#F090A0", "#50D050",\n "#C88033", "#7D80B0", "#C28F8F", "#668F8F", "#BD80E3", "#FFA100", "#A62929",\n "#5CB8D1", "#702EB0", "#00FF00", "#94FFFF", "#94E0E0", "#73C2C9", "#54B5B5",\n "#3B9E9E", "#248F8F", "#0A7D8C", "#006985", "#C0C0C0", "#FFD98F", "#A67573",\n "#668080", "#9E63B5", "#D47A00", "#940094", "#429EB0", "#57178F", "#00C900",\n "#70D4FF", "#FFFFC7", "#D9FFC7", "#C7FFC7", "#A3FFC7", "#8FFFC7", "#61FFC7",\n "#45FFC7", "#30FFC7", "#1FFFC7", "#00FF9C", "#00E675", "#00D452", "#00BF38",\n "#00AB24", "#4DC2FF", "#4DA6FF", "#2194D6", "#267DAB", "#266696", "#175487",\n "#D0D0E0", "#FFD123", "#B8B8D0", "#A6544D", "#575961", "#9E4FB5", "#AB5C00",\n "#754F45", "#428296", "#420066", "#007D00", "#70ABFA", "#00BAFF", "#00A1FF",\n "#008FFF", "#0080FF", "#006BFF", "#545CF2", "#785CE3", "#8A4FE3", "#A136D4",\n "#B31FD4", "#B31FBA", "#B30DA6", "#BD0D87", "#C70066", "#CC0059", "#D1004F",\n "#D90045", "#E00038", "#E6002E", "#EB0026")\nnames(element_colors) <- c("H", "He", "Li", "Be", "B", "C", "N", "O", "F", "Ne", "Na", "Mg", "Al", "Si",\n "P", "S", "Cl", "Ar", "K", "Ca", "Sc", "Ti", "V", "Cr", "Mn", "Fe", "Co", "Ni",\n "Cu", "Zn", "Ga", "Ge", "As", "Se", "Br", "Kr", "Rb", "Sr", "Y", "Zr", "Nb", "Mo",\n "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", "In", "Sn", "Sb", "Te", "I", "Xe", "Cs", "Ba",\n "La", "Ce", "Pr", "Nd", "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb",\n "Lu", "Hf", "Ta", "W", "Re", "Os", "Ir", "Pt", "Au", "Hg", "Tl", "Pb", "Bi", "Po",\n "At", "Rn", "Fr", "Ra", "Ac", "Th", "Pa", "U", "Np", "Pu", "Am", "Cm", "Bk", "Cf",\n "Es", "Fm", "Md", "No", "Lr", "Rf", "Db", "Sg", "Bh", "Hs", "Mt")\n\nxray_frame <- as.data.frame(xray_brick, xy=TRUE) \nxray_frame <- gather(xray_frame, element, value, colnames(xray_frame)[3]:colnames(xray_frame)[ncol(xray_frame)])\n\n\n\nelement_plotter<-function(coord_frame, brick, SEM_image, colors){\n p <-rasterVis::gplot(SEM_image) +\n geom_tile(aes(fill = value)) +\n scale_fill_gradient(low = 'black', high = 'white') +\n ggnewscale::new_scale_fill()\n for(i in names(brick)){ \n message(paste0("Adding ", names(brick[[i]]), " to plot..."))\n element_coords <- coord_frame %>%\n filter(element == names(brick[[i]]) & value!=0)\n p <- p+geom_raster(element_coords, mapping = aes(x, y, fill = element, alpha = value)) +\n scale_fill_manual(values = colors) +\n ggnewscale::new_scale_fill()\n }\n message("Writing plot...")\n suppressWarnings(print(p + \n coord_fixed() +\n theme(axis.title = element_blank(),\n axis.text = element_blank(),\n legend.position = "none")))\n}\n\nelement_plot <- element_plotter(xray_frame, xray_brick, SEM_panorama_merged, element_colors)\n\nelement_plot_legend <- data.frame(element = unique(xray_frame$element)) %>%\n rownames_to_column() %>% \n ggplot(aes(element, rowname, fill=element)) + \n geom_bar(stat= "identity") + \n scale_fill_manual(values = element_colors) +\n guides(fill=guide_legend(nrow=2)) +\n theme(legend.title = ggplot2::element_blank(),\n legend.text = ggplot2::element_text(size = 8))\n\nelement_plot_with_legend <- plot_grid(\n element_plot, \n plot_grid(get_legend(element_plot_legend), \n ncol = 1), \n nrow = 2, \n rel_heights = c(8,2)\n)\n\n\n# generate PDF ------------------------------------------------------------\n\nif(!is.na(opt$n)){\n if(!is.na(opt$o)){\n pdf(paste0(opt$o, "/", opt$n, "_element_plot.pdf"),\n width = 13.33, \n height = 7.5)\n }else{\n pdf(paste0(opt$n, "_element_plot.pdf"),\n width = 13.33, \n height = 7.5)\n }\n}else{\n if(!is.na(opt$o)){\n pdf(paste0(opt$o, "/", "element_plot.pdf"),\n width = 13.33, \n height = 7.5)\n }else{\n pdf("element_plot.pdf",\n width = 13.33, \n height = 7.5) \n }\n}\n\n\nprint(element_plot_with_legend)\n\ndev.off()\n\nmessage("...complete.")\n\nremove(list = ls())\n\n}\n\n