diff --git a/NEWS.md b/NEWS.md index 036500fa..2fdeef9c 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,13 +1,15 @@ # VoltRon 0.2.7 -- `transferData` now allows integrating single cell data object (`Seurat` or `SingleCellExperiment`) - for transfering features (e.g. gene expression) or metadata features (e.g. cell types, annotations). - `registerSpatialData` now allows registering assays with no images. In case of either one of the assays, do not have images, assays can be registered with only the manual approach. +- The Shiny interface for `registerSpatialData` now reports on the accuracy of the alignment using multiple + metrics such as Intersection, Bhattacharyya and Matte's Mutual Information. - An image-free alignment tutorial has been added where DBIT-Seq and a QuPath processed mIF experiment are aligned using manually selected landmarks. +- `transferData` now allows integrating single cell data object (`Seurat` or `SingleCellExperiment`) + for transfering features (e.g. gene expression) or metadata features (e.g. cell types, annotations). - Now `importImageData` and `importQuPathIF` functions only work with segments already converted to a - list by the user, or the `generateSegments` function whose arguement is an **sf** object. + list by the user, or the `generateSegments` function whose argument is an **sf** object. - The `formAssay` function now accepts only segments without the user manually generating coordinates (or centroids). - The `name` argument in functions like `vrImages`, `vrCoordinates` and `vrSegments` is replaced diff --git a/R/RcppExports.R b/R/RcppExports.R index 8ab45350..e534540a 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -1,8 +1,12 @@ # Generated by using Rcpp::compileAttributes() -> do not edit by hand # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 -automated_registeration_rawvector <- function(ref_image, query_image, width1, height1, width2, height2, GOOD_MATCH_PERCENT, MAX_FEATURES, invert_query, invert_ref, flipflop_query, flipflop_ref, rotate_query, rotate_ref, matcher, method, nonrigid) { - .Call('_VoltRon_automated_registeration_rawvector', PACKAGE = 'VoltRon', ref_image, query_image, width1, height1, width2, height2, GOOD_MATCH_PERCENT, MAX_FEATURES, invert_query, invert_ref, flipflop_query, flipflop_ref, rotate_query, rotate_ref, matcher, method, nonrigid) +accuracy_rawvector <- function(ref_image, query_image, mask, width, height, type, overlay_images = TRUE, compute_matte_map = TRUE) { + .Call('_VoltRon_accuracy_rawvector', PACKAGE = 'VoltRon', ref_image, query_image, mask, width, height, type, overlay_images, compute_matte_map) +} + +automated_registeration_rawvector <- function(ref_image, query_image, width1, height1, width2, height2, GOOD_MATCH_PERCENT, MAX_FEATURES, invert_query, invert_ref, flipflop_query, flipflop_ref, rotate_query, rotate_ref, matcher, method, nonrigid, compute_matte_map = TRUE) { + .Call('_VoltRon_automated_registeration_rawvector', PACKAGE = 'VoltRon', ref_image, query_image, width1, height1, width2, height2, GOOD_MATCH_PERCENT, MAX_FEATURES, invert_query, invert_ref, flipflop_query, flipflop_ref, rotate_query, rotate_ref, matcher, method, nonrigid, compute_matte_map) } replaceNaMatrix <- function(mat, replace) { @@ -21,8 +25,8 @@ warpImageManual <- function(ref_image, query_image, mapping, width1, height1, wi .Call('_VoltRon_warpImageManual', PACKAGE = 'VoltRon', ref_image, query_image, mapping, width1, height1, width2, height2) } -manual_registeration_rawvector <- function(ref_image, query_image, reference_landmark, query_landmark, width1, height1, width2, height2, method, nonrigid) { - .Call('_VoltRon_manual_registeration_rawvector', PACKAGE = 'VoltRon', ref_image, query_image, reference_landmark, query_landmark, width1, height1, width2, height2, method, nonrigid) +manual_registeration_rawvector <- function(ref_image, query_image, reference_landmark, query_landmark, width1, height1, width2, height2, invert_query, invert_ref, method, nonrigid, compute_matte_map = TRUE) { + .Call('_VoltRon_manual_registeration_rawvector', PACKAGE = 'VoltRon', ref_image, query_image, reference_landmark, query_landmark, width1, height1, width2, height2, invert_query, invert_ref, method, nonrigid, compute_matte_map) } manual_registeration_matrix <- function(query_data, reference_landmark, query_landmark, method, nonrigid) { @@ -33,6 +37,10 @@ applyRcppMapping <- function(coords, mapping) { .Call('_VoltRon_applyRcppMapping', PACKAGE = 'VoltRon', coords, mapping) } +generateOverlapMask <- function(dsize, trans_mat, ssize) { + .Call('_VoltRon_generateOverlapMask', PACKAGE = 'VoltRon', dsize, trans_mat, ssize) +} + build_snn_rank <- function(neighbors) { .Call('_VoltRon_build_snn_rank', PACKAGE = 'VoltRon', neighbors) } diff --git a/R/auxiliary.R b/R/auxiliary.R index 71ee0a62..7af6232a 100644 --- a/R/auxiliary.R +++ b/R/auxiliary.R @@ -58,6 +58,43 @@ fixVoltRon <- function(object) { object } +#### +# Bioformats extensions for image pyramids #### +#### + +.PYRAMID_FORMATS <- c( + ".ome.tiff", # OME-TIFF pyramidal TIFF variants (OME-TIFF and related) + ".ome.tif", + ".ome.tf2", + ".ome.tf8", + ".ome.btf", + ".svs", # Aperio SVS + ".afi", # Aperio AFI + ".ndpi", # Hamamatsu NDPI + ".ndpis", # Hamamatsu NDPI variants + ".qptiff" # CODEX ? +) + +#### +# Alignment Metrics #### +#### + +.ALIGNMENT_ACCURACY_METRICS <- c( + "Intersection", + "Bhattacharyya", + "Matte's MI" +) + +.ALIGNMENT_KEYPOINT_METRICS <- c( + "#Keypoints", + "Inlier Perc.", + "sd ref. kpts (>1?)", + "sd query kpts (>1?)", + "sd grid (in [w,h]?)", + "Median distance", + "Degenerate" +) + #### # Matrix Operations #### #### @@ -326,24 +363,6 @@ getBasilisk <- function() { py_env } -#### -# Bioformats extensions for pyramids #### -#### - -.PYRAMID_FORMATS <- c( - ".ome.tiff", # OME-TIFF pyramidal TIFF variants (OME-TIFF and related) - ".ome.tif", - ".ome.tf2", - ".ome.tf8", - ".ome.btf", - ".svs", # Aperio SVS - ".afi", # Aperio AFI - ".ndpi", # Hamamatsu NDPI - ".ndpis", # Hamamatsu NDPI variants - ".qptiff" # CODEX ? -) - - #### # Other Auxiliary tools #### #### diff --git a/R/registration.R b/R/registration.R index 1f5e9a52..c56ce166 100644 --- a/R/registration.R +++ b/R/registration.R @@ -619,12 +619,23 @@ getAlignmentTabPanel <- function(len_images, centre, register_ind) { do.call( tabsetPanel, c( + # alignment tab panel id = 'image_tab_panel_alignment', lapply(register_ind, function(i) { tabPanel( paste0("Ali. ", i, "->", centre), br(), - fluidRow(imageOutput(paste0("plot_alignment", i))) + + tabsetPanel( + id = paste0("inner_tabs", i), + tabPanel("Alignment Stat.", + tableOutput(paste0("alignment_stats", i))), + tabPanel("Matte's MI Map", + imageOutput(paste0("plot_matte_map", i))), + tabPanel("Matching Keypoints", + imageOutput(paste0("plot_keypoint_match", i))) + + ) ) }) ) @@ -833,7 +844,7 @@ updateParameterPanels <- function(len_images, params, input, output, session) { "nonrigid", choices = c( "TPS (OpenCV)", - "BSpline (SimpleITK)" + if(input$Method != "Non-Rigid") "BSpline (SimpleITK)" else NULL ), selected = "TPS (OpenCV)" ) @@ -1387,94 +1398,6 @@ applyPerspectiveTransform <- function( return(object) } -#### -# Managing Mappings #### -#### - -manageMapping <- function(mappings) { - # check if all transformations are homography - allHomography <- suppressWarnings(all(lapply(mappings, function(map) { - nrow(map[[1]] > 0) && is.null(map[[2]]) - }))) - - # change the mapping - new_mappings <- list() - if (allHomography) { - mappings <- lapply(mappings, function(map) map[[1]]) - new_mappings <- list( - list(Reduce("%*%", mappings), NULL) - ) - } else { - new_mappings <- mappings - } - - # return - return(new_mappings) -} - -applyMapping <- function(coords, mapping){ - mapping_new <- mapping - if(!is.null(mapping[[1]][[2]])){ - if(is(mapping[[1]][[2]][[1]], "_p_itk__simple__TransformixImageFilter")){ - mapping_new[[1]] <- list(mapping[[1]][[1]], NULL) - coords <- applyRcppMapping(coords, mapping_new) - coords <- applySimpleITKMapping(coords, mapping[[1]][[2]][[1]]) - } else { - coords <- applyRcppMapping(coords, mapping) - } - } else { - coords <- applyRcppMapping(coords, mapping) - } - coords -} - -#' @importFrom utils write.table -applySimpleITKMapping <- function(coords, mapping){ - - # check SimpleITK - if (!requireNamespace('SimpleITK')) { - stop("Please install SimpleITK package!: ", - "remotes::install_github('BIMSBbioinfo/SimpleITKRInstaller')", - ", this is gonna take a while :)") - } - - # temp dir, delete later - tmpdir <- tempdir() - tmpdir <- file.path(tmpdir, "SimpleITK") - dir.create(tmpdir, showWarnings = FALSE) - - # get image - input_file <- file.path(tmpdir, "inputpoints.txt") - output_file <- file.path(tmpdir, "outputpoints.txt") - - # apply transformation - tfx <- mapping - suppressWarnings(file.remove(input_file, showWarnings = FALSE)) - cat("point\n", nrow(coords), "\n", file = input_file) - utils::write.table(coords, input_file, append = TRUE, - col.names = FALSE, row.names = FALSE, quote = FALSE) - tfx$SetOutputDirectory(tmpdir) - tfx$SetFixedPointSetFileName(input_file) - tmp <- tfx$Execute() - - # get points - lines_coords <- readLines(output_file) - coords <- do.call( - rbind, - lapply(lines_coords, function(x) { - tmp <- strsplit(strsplit(x, split = "\\t")[[1]][6], - split = " ")[[1]][c(5,6)] - as.numeric(tmp) - }) - ) - - # delete dir - unlink(tmpdir, recursive = TRUE) - - # return - coords -} - #### # Managing Parameters #### #### @@ -2555,6 +2478,111 @@ ggplot_to_magick <- function(plot, extent = NULL, width = 8, height = 6, dpi = 3 magick::image_read(tf) } +#' @noRd +convertToSitkImage <- function(img){ + # img_info <- magick::image_info(img) + img_data <- magick::image_data(img, channels = "gray") + dim_img <- dim(img_data) + img_data <- as.vector(img_data, mode = "integer") + # dim(img_data) <- c(3, img_info$width, img_info$height) + dim(img_data) <- dim_img + img_data <- aperm(img_data, c(2,3,1)) + SimpleITK::as.image(img_data, isVector = TRUE) +} + +#### +# Managing Mappings #### +#### + +manageMapping <- function(mappings) { + # check if all transformations are homography + allHomography <- suppressWarnings(all(lapply(mappings, function(map) { + nrow(map[[1]] > 0) && is.null(map[[2]]) + }))) + + # change the mapping + new_mappings <- list() + if (allHomography) { + mappings <- lapply(mappings, function(map) map[[1]]) + new_mappings <- list( + list(Reduce("%*%", mappings), NULL) + ) + } else { + new_mappings <- mappings + } + + # return + return(new_mappings) +} + +applyMapping <- function(coords, mapping){ + mapping_new <- mapping + if(!is.null(mapping[[1]][[2]])){ + if(is(mapping[[1]][[2]][[1]], "_p_itk__simple__TransformixImageFilter")){ + mapping_new[[1]] <- list(mapping[[1]][[1]], NULL) + coords <- applyRcppMapping(coords, mapping_new) + coords <- applySimpleITKMapping(coords, mapping[[1]][[2]][[1]]) + } else { + coords <- applyRcppMapping(coords, mapping) + } + } else { + coords <- applyRcppMapping(coords, mapping) + } + coords +} + +#' @importFrom utils write.table +applySimpleITKMapping <- function(coords, mapping){ + + # check SimpleITK + if (!requireNamespace('SimpleITK')) { + stop("Please install SimpleITK package!: ", + "remotes::install_github('BIMSBbioinfo/SimpleITKRInstaller')", + ", this is gonna take a while :)") + } + + # temp dir, delete later + tmpdir <- tempdir() + tmpdir <- file.path(tmpdir, "SimpleITK") + dir.create(tmpdir, showWarnings = FALSE) + + # get image + input_file <- file.path(tmpdir, "inputpoints.txt") + output_file <- file.path(tmpdir, "outputpoints.txt") + + # apply transformation + tfx <- mapping + suppressWarnings(file.remove(input_file, showWarnings = FALSE)) + cat("point\n", nrow(coords), "\n", file = input_file) + utils::write.table(coords, input_file, append = TRUE, + col.names = FALSE, row.names = FALSE, quote = FALSE) + tfx$SetOutputDirectory(tmpdir) + tfx$SetFixedPointSetFileName(input_file) + tmp <- tfx$Execute() + + # get points + lines_coords <- readLines(output_file) + coords <- do.call( + rbind, + lapply(lines_coords, function(x) { + tmp <- strsplit(strsplit(x, split = "\\t")[[1]][6], + split = " ")[[1]][c(5,6)] + as.numeric(tmp) + }) + ) + + # delete dir + unlink(tmpdir, recursive = TRUE) + + # return + coords +} + + +#### +# Managing Transformations #### +#### + #' transformImage #' #' Apply given transformations to a magick image @@ -2697,7 +2725,6 @@ warpImage <- function(ref_image, query_image, mapping) { if(is(mapping[[1]][[2]][[2]], "_p_itk__simple__TransformixImageFilter")){ query_image <- magick::image_read(query_image) query_image <- warpSimpleITKImage( - ref_image = ref_image, query_image = query_image, mapping = mapping[[1]][[2]][[2]] ) @@ -2737,14 +2764,13 @@ warpImage <- function(ref_image, query_image, mapping) { #' #' Warping a query image given a homography image #' -#' @param ref_image reference image #' @param query_image query image #' @param mapping a list of the homography matrices and TPS keypoints #' #' @importFrom magick image_read image_data #' #' @export -warpSimpleITKImage <- function(ref_image, query_image, mapping) { +warpSimpleITKImage <- function(query_image, mapping) { # check SimpleITK if (!requireNamespace('SimpleITK')) { @@ -2862,6 +2888,8 @@ getManualRegisteration <- function( { # Register keypoints aligned_image_list <- list() + matte_map_list <- list() + alignment_stats_list <- list() for (i in register_ind) { # Increment the progress bar, and update the detail text. incProgress( @@ -2883,6 +2911,12 @@ getManualRegisteration <- function( # save matches aligned_image_list[[i]] <- results$aligned_image + + # save matte map + matte_map_list[[i]] <- results$matte_map + + # save alignment stats + alignment_stats_list[[i]] <- results$alignment_stats } } ) @@ -2930,7 +2964,45 @@ getManualRegisteration <- function( deleteFile = TRUE ) }) + + # Plot Matte + lapply(register_ind, function(i) { + if(length(matte_map_list)){ + cur_alignment_image <- matte_map_list[[i]] + output[[paste0("plot_matte_map", i)]] <- renderPlot({ + if (!suppressWarnings(!is.matrix(cur_alignment_image))) { + cur_alignment_image <- + cur_alignment_image[nrow(cur_alignment_image):1,] + cur_alignment_image <- + as.data.frame(as.table(cur_alignment_image)) + ggplot(cur_alignment_image, + aes(Var2, Var1, fill= Freq)) + + ggplot2::geom_tile() + + ggplot2::theme_void() + + ggplot2::coord_fixed(expand = FALSE) + + ggplot2::scale_fill_gradient(low = "#440154FF", + high = "#FDE725FF", + # na.value = NA, + name = "Matte's MI") + } + }) + } + }) + # Plot Alignment Stats + lapply(register_ind, function(i) { + if(length(alignment_stats_list)){ + cur_align_stats <- alignment_stats_list[[i]] + output[[paste0("alignment_stats", i)]] <- renderTable({ + tab <- data.frame(Metrics = names(cur_align_stats[["coarse"]]), + `Coarse` = cur_align_stats[["coarse"]]) + if(!all(is.na(cur_align_stats[["fine"]]))) + tab$Fine <- cur_align_stats[["fine"]] + tab + }, digits = 5, na = "") + } + }) + # Output summary output[["summary"]] <- renderUI({ str1 <- paste0(" Registration Summary:") @@ -2959,7 +3031,8 @@ computeManualPairwiseTransform <- function( keypoints_list, query_ind, ref_ind, - input + input, + compute_matte_map = TRUE ) { # determine the number of transformation to map from query to the reference indices <- query_ind:ref_ind @@ -2998,8 +3071,23 @@ computeManualPairwiseTransform <- function( ref_image, target_landmark, reference_landmark, + invert_query = input[[paste0( + "negate_", + query_label, + "_image", + cur_map[1] + )]] == + "Yes", + invert_ref = input[[paste0( + "negate_", + ref_label, + "_image", + cur_map[2] + )]] == + "Yes", method = input$Method, - nonrigid = if(is.null(input$nonrigid)) "None" else input$nonrigid + nonrigid = if(is.null(input$nonrigid)) "None" else input$nonrigid, + compute_matte_map = compute_matte_map ) # run SimpleITK as fine registration @@ -3015,25 +3103,45 @@ computeManualPairwiseTransform <- function( tfx <- getSimpleITKAutomatedRegistration( ref_image = ref_image, query_image = query_image, - invert_query = FALSE, - invert_ref = FALSE, + invert_query = input[[paste0( + "negate_", + query_label, + "_image", + cur_map[1] + )]] == + "Yes", + invert_ref = input[[paste0( + "negate_", + ref_label, + "_image", + cur_map[2] + )]] == + "Yes", flipflop_query = FALSE, flipflop_ref = FALSE, rotate_query = FALSE, rotate_ref = FALSE, - initial_mapping = list(reg[[1]]) + initial_mapping = list(reg[[1]]), + compute_matte_map = compute_matte_map ) - reg$aligned_image <- tfx$aligned_image reg[[1]][[2]] <- tfx$transformation + reg$aligned_image <- tfx$aligned_image + reg$matte_map <- tfx$matte_map + reg$alignment_stats$fine[names(tfx$alignment_metrics)] <- + tfx$alignment_metrics } # return transformation matrix and images mapping[[kk]] <- reg[[1]] aligned_image <- reg$aligned_image + matte_map <- reg$matte_map + alignment_stats <- reg$alignment_stats } return(list(mapping = mapping, - aligned_image = aligned_image)) + aligned_image = aligned_image, + matte_map = matte_map, + alignment_stats = alignment_stats)) } #' getRcppManualRegistration @@ -3055,6 +3163,8 @@ getRcppManualRegistration <- function( ref_image, query_landmark, reference_landmark, + invert_query = FALSE, + invert_ref = FALSE, method = "Homography", nonrigid = "TPS (OpenCV)" ) { @@ -3089,7 +3199,6 @@ getRcppManualRegistration <- function( query_landmark[, 2] <- dim(query_image)[3] - query_landmark[, 2] } - reg <- if(ncol(query_image) == 2){ manual_registeration_matrix( @@ -3105,6 +3214,8 @@ getRcppManualRegistration <- function( query_image, reference_landmark = reference_landmark, query_landmark = query_landmark, + invert_query = invert_query, + invert_ref = invert_ref, width1 = dim(ref_image)[2], height1 = dim(ref_image)[3], width2 = dim(query_image)[2], @@ -3133,9 +3244,55 @@ getRcppManualRegistration <- function( magick::image_read(reg[[2]]) } + # check for null data + if(length(reg) > 2){ + + # check matte map + matte_map <- + if (!is.null(reg[[3]])) { + tmp <- reg[[3]] + tmp[tmp < 0] <- 0 + tmp + } else NA + + # check alignment statistics + alignment_stats <- list() + metrics <- .ALIGNMENT_ACCURACY_METRICS + metrics_set <- setNames(rep(NA, length(metrics)), metrics) + alignment_stats[["coarse"]] <- { + if (!is.null(reg[[4]])){ + if(!all(names(reg[[4]]) %in% metrics)){ + stop("There are missing accuracy metrics!") + } else { + metrics_set[metrics] <- reg[[4]][metrics] + metrics_set + } + } else{ + NA + } + } + alignment_stats[["fine"]] <- { + if (!is.null(reg[[5]])){ + if(!all(names(reg[[5]]) %in% metrics)){ + stop("There are missing accuracy metrics!") + } else { + metrics_set[metrics] <- reg[[5]][metrics] + metrics_set + } + } else{ + NA + } + } + } else { + matte_map <- NULL + alignment_stats <- NULL + } + return(list( transmat = reg[[1]], - aligned_image = aligned_image + aligned_image = aligned_image, + matte_map = matte_map, + alignment_stats = alignment_stats )) } @@ -3195,6 +3352,8 @@ getAutomatedRegisteration <- function( overlayed_image_list <- list() aligned_image_list <- list() alignment_image_list <- list() + matte_map_list <- list() + alignment_stats_list <- list() for (i in register_ind) { # Increment the progress bar, and update the detail text. incProgress( @@ -3225,6 +3384,12 @@ getAutomatedRegisteration <- function( # save matches alignment_image_list[[i]] <- results$alignment_image + + # save matte map + matte_map_list[[i]] <- results$matte_map + + # save alignment stats + alignment_stats_list[[i]] <- results$alignment_stats } } ) @@ -3257,12 +3422,46 @@ getAutomatedRegisteration <- function( # Plot Alignment lapply(register_ind, function(i) { cur_alignment_image <- alignment_image_list[[i]] - output[[paste0("plot_alignment", i)]] <- renderPlot({ + output[[paste0("plot_keypoint_match", i)]] <- renderPlot({ if (!suppressWarnings(is.na(cur_alignment_image))) { magick::image_ggplot(cur_alignment_image) } }) }) + + # Plot Matte + lapply(register_ind, function(i) { + cur_alignment_image <- matte_map_list[[i]] + output[[paste0("plot_matte_map", i)]] <- renderPlot({ + if (!suppressWarnings(!is.matrix(cur_alignment_image))) { + cur_alignment_image <- + cur_alignment_image[nrow(cur_alignment_image):1,] + cur_alignment_image <- + as.data.frame(as.table(cur_alignment_image)) + ggplot(cur_alignment_image, + aes(Var2, Var1, fill= Freq)) + + ggplot2::geom_tile() + + ggplot2::theme_void() + + ggplot2::coord_fixed(expand = FALSE) + + ggplot2::scale_fill_gradient(low = "#440154FF", + high = "#FDE725FF", + # na.value = NA, + name = "Matte's MI") + } + }) + }) + + # Plot Alignment Stats + lapply(register_ind, function(i) { + cur_align_stats <- alignment_stats_list[[i]] + output[[paste0("alignment_stats", i)]] <- renderTable({ + tab <- data.frame(Metrics = names(cur_align_stats[["coarse"]]), + `Coarse` = cur_align_stats[["coarse"]]) + if(!all(is.na(cur_align_stats[["fine"]]))) + tab$Fine <- cur_align_stats[["fine"]] + tab + }, digits = 5, na = "") + }) # Output summary output[["summary"]] <- renderUI({ @@ -3292,7 +3491,8 @@ computeAutomatedPairwiseTransform <- function( channel_names, query_ind, ref_ind, - input + input, + compute_matte_map = TRUE ) { # determine the number of transformation to map from query to the reference indices <- query_ind:ref_ind @@ -3384,7 +3584,8 @@ computeAutomatedPairwiseTransform <- function( rotate_ref = input[[paste0("rotate_", ref_label, "_image", cur_map[2])]], matcher = input$Matcher, method = input$Method, - nonrigid = if(is.null(input$nonrigid)) "None" else input$nonrigid + nonrigid = if(is.null(input$nonrigid)) "None" else input$nonrigid, + compute_matte_map = compute_matte_map ) # update transformation matrix @@ -3451,10 +3652,15 @@ computeAutomatedPairwiseTransform <- function( )]], rotate_ref = input[[paste0( "rotate_", ref_label, "_image", cur_map[2])]], - initial_mapping = list(reg[[1]]) + initial_mapping = list(reg[[1]]), + compute_matte_map = compute_matte_map ) - reg$aligned_image <- tfx$aligned_image reg[[1]][[2]] <- tfx$transformation + reg$aligned_image <- tfx$aligned_image + reg$overlay_image <- tfx$overlay_image + reg$matte_map <- tfx$matte_map + reg$alignment_stats$fine[names(tfx$alignment_metrics)] <- + tfx$alignment_metrics } # return transformation matrix and images @@ -3463,6 +3669,8 @@ computeAutomatedPairwiseTransform <- function( aligned_image <- reg$aligned_image alignment_image <- reg$alignment_image overlay_image <- reg$overlay_image + matte_map <- reg$matte_map + alignment_stats <- reg$alignment_stats } return(list( @@ -3470,7 +3678,9 @@ computeAutomatedPairwiseTransform <- function( dest_image = dest_image, aligned_image = aligned_image, alignment_image = alignment_image, - overlay_image = overlay_image + overlay_image = overlay_image, + matte_map = matte_map, + alignment_stats = alignment_stats )) } @@ -3508,7 +3718,8 @@ getRcppAutomatedRegistration <- function( rotate_ref = "0", matcher = "FLANN", method = "Homography", - nonrigid = "TPS (OpenCV)" + nonrigid = "TPS (OpenCV)", + compute_matte_map = TRUE ) { ref_image <- magick::image_data(ref_image, channels = "rgb") query_image <- magick::image_data(query_image, channels = "rgb") @@ -3530,7 +3741,8 @@ getRcppAutomatedRegistration <- function( rotate_ref = rotate_ref, matcher = matcher, method = method, - nonrigid = nonrigid + nonrigid = nonrigid, + compute_matte_map = compute_matte_map ) # check for null keypoints @@ -3538,13 +3750,52 @@ getRcppAutomatedRegistration <- function( reg[[1]] <- list(reg[[1]][[1]], NULL) } - # check for failed registeration + # check for failed registration aligned_image <- if (!is.null(reg[[3]])) magick::image_read(reg[[3]]) else NA alignment_image <- if (!is.null(reg[[4]])) magick::image_read(reg[[4]]) else NA overlay_image <- if (!is.null(reg[[5]])) magick::image_read(reg[[5]]) else NA + + # check matte maps + matte_map <- + if (!is.null(reg[[6]])){ + tmp <- reg[[6]] + tmp[tmp < 0] <- 0 + tmp + + } else NA + + # check alignment statistics + alignment_stats <- list() + metrics <- c(.ALIGNMENT_ACCURACY_METRICS, + .ALIGNMENT_KEYPOINT_METRICS) + metrics_set <- setNames(rep(NA, length(metrics)), metrics) + alignment_stats[["coarse"]] <- { + if (!is.null(reg[[7]])){ + if(!all(names(reg[[7]]) %in% metrics)){ + stop("There are missing accuracy metrics!") + } else { + metrics_set[metrics] <- reg[[7]][metrics] + metrics_set + } + } else{ + NA + } + } + alignment_stats[["fine"]] <- { + if (!is.null(reg[[8]])){ + if(!all(names(reg[[8]]) %in% metrics)){ + stop("There are missing accuracy metrics!") + } else { + metrics_set[metrics] <- reg[[8]][metrics] + metrics_set + } + } else{ + NA + } + } # return return(list( @@ -3552,7 +3803,9 @@ getRcppAutomatedRegistration <- function( dest_image = magick::image_read(reg[[2]]), aligned_image = aligned_image, alignment_image = alignment_image, - overlay_image = overlay_image + overlay_image = overlay_image, + matte_map = matte_map, + alignment_stats = alignment_stats )) } @@ -3582,7 +3835,8 @@ getSimpleITKAutomatedRegistration <- function( flipflop_ref = "None", rotate_query = "0", rotate_ref = "0", - initial_mapping = NULL + initial_mapping = NULL, + compute_matte_map = TRUE ){ # check SimpleITK if (!requireNamespace('SimpleITK')) { @@ -3613,7 +3867,6 @@ getSimpleITKAutomatedRegistration <- function( ref_image <- array(as.raw(ref_image), dim = dim(ref_image)) ref_image <- magick::image_read(ref_image) } - query_image <- rotateImage(query_image, as.numeric(rotate_query)) if (flipflop_query == "Flip") { query_image <- flipImage(query_image) @@ -3622,39 +3875,44 @@ getSimpleITKAutomatedRegistration <- function( } if(invert_query) query_image <- negateImage(query_image) - # query_image <- magick::image_negate(query_image) + + # generate coarse mapped mask by warping + ref_info <- getImageInfo(ref_image) + query_info <- getImageInfo(query_image) + mask <- generateOverlapMask(c(ref_info$width, ref_info$height), + initial_mapping[[1]][[1]], + c(query_info$width, query_info$height)) + # mask <- magick::image_read(mask) + + # warp image query_image <- warpImage(ref_image = ref_image, query_image = query_image, mapping = initial_mapping) - # prepare images - ref_image1 <- magick::as_EBImage(ref_image) - # ref_image1 <- EBImage::imageData(ref_image1) - # dim_img <- 1:length(dim(ref_image)) - # dim_img[1:2] <- rev(dim_img[1:2]) - # ref_image1 <- aperm(ref_image1, perm = c(2,1,3)) - EBImage::writeImage(ref_image1, - files = file.path(tmpdir, "ref_image.tiff"), - compression = "LZW", reduce = TRUE) - fixed <- SimpleITK::ReadImage(file.path(tmpdir, "ref_image.tiff"), - 'sitkUInt8') - query_image1 <- as_EBImage(query_image) - # dim_img <- 1:length(dim(query_image)) - # dim_img[1:2] <- rev(dim_img[1:2]) - # query_image1 <- EBImage::imageData(query_image1) - # query_image1 <- aperm(query_image1, perm = c(2,1)) - EBImage::writeImage(query_image1, - files = file.path(tmpdir, "query_image.tiff"), - compression = "LZW", reduce = TRUE) - moving <- SimpleITK::ReadImage(file.path(tmpdir, "query_image.tiff"), - 'sitkUInt8') + # prepare images and masks + # magick::image_write(ref_image, file.path(tmpdir, "ref_image.tiff"), + # compression = "LZW") + # fixed <- SimpleITK::ReadImage(file.path(tmpdir, "ref_image.tiff"), + # 'sitkUInt8') + # magick::image_write(query_image, file.path(tmpdir, "query_image.tiff"), + # compression = "LZW") + # moving <- SimpleITK::ReadImage(file.path(tmpdir, "query_image.tiff"), + # 'sitkUInt8') + fixed <- convertToSitkImage(ref_image) + fixed <- SimpleITK::Cast(fixed, "sitkUInt8") + moving <- convertToSitkImage(query_image) + moving <- SimpleITK::Cast(moving, "sitkUInt8") + mask <- SimpleITK::as.image(array(mask, rev(dim(mask)))) + # mask <- SimpleITK::as.image(array(as.integer(mask != 0L), + # rev(dim(mask)))) + mask <- SimpleITK::Cast(mask, "sitkUInt8") # get registration for image elx <- SimpleITK::ElastixImageFilter() elx$SetOutputDirectory(tmpdir) elx$SetFixedImage(fixed) elx$SetMovingImage(moving) - parameterMapVector = SimpleITK::VectorOfParameterMap() + # elx$SetMovingMask(mask) mp <- SimpleITK:::ReadParameterFile( system.file("extdata", "bspline_map.txt", package = "VoltRon") ) @@ -3663,21 +3921,33 @@ getSimpleITKAutomatedRegistration <- function( tmp <- elx$Execute() sitk_img <- SimpleITK::ReadImage(file.path(tmpdir, "result.0.tif")) arr <- SimpleITK::as.array(sitk_img) - arr8 <- 255 * (arr - min(arr)) / (max(arr) - min(arr)) - arr8 <- array(as.integer(arr8), dim = dim(arr)) - arr8 <- aperm(arr8, perm = c(2,1)) - aligned_image <- magick::image_read(as.raster(arr8 / 255)) + arr <- (arr - min(arr)) / (max(arr) - min(arr)) + arr <- array(arr, dim = dim(arr)) + arr <- aperm(arr, perm = c(2,1)) + aligned_image <- magick::image_read(as.raster(arr)) transform_param_map <- elx$GetTransformParameterMap() tfx_image <- SimpleITK::TransformixImageFilter() tfx_image$LogToConsoleOff() tfx_image$SetTransformParameterMap(transform_param_map) + + # warp mask + tfx_mask <- SimpleITK::TransformixImageFilter() + tfx_mask$LogToConsoleOff() + tfx_mask$SetTransformParameterMap(transform_param_map) + mask$CopyInformation(moving) + tfx_mask$SetMovingImage(mask) + tmp <- tfx_mask$Execute() + aligned_mask <- SimpleITK::as.array(tfx_mask$GetResultImage()) + aligned_mask <- array(aligned_mask, dim = c(dim(aligned_mask), 1)) + aligned_mask <- aperm(aligned_mask, c(2,1,3)) + aligned_mask <- magick::image_read(aligned_mask) # get transformation for the points and observations elx <- SimpleITK::ElastixImageFilter() elx$SetOutputDirectory(tmpdir) elx$SetFixedImage(moving) + elx$SetFixedMask(mask) elx$SetMovingImage(fixed) - parameterMapVector = SimpleITK::VectorOfParameterMap() mp <- SimpleITK:::ReadParameterFile( system.file("extdata", "bspline_map.txt", package = "VoltRon") ) @@ -3693,8 +3963,30 @@ getSimpleITKAutomatedRegistration <- function( # delete dir unlink(tmpdir, recursive = TRUE) + # calculate alignment accuracy + results <- getAlignmentAccuracy(ref_image, + aligned_image, + aligned_mask, + "Fine", + compute_matte_map) + + # convert images + overlay_image <- + if (!is.null(results[[3]])) magick::image_read(results[[3]]) else NA + + # check matte maps + matte_map <- + if (!is.null(results[[2]])){ + tmp <- results[[2]] + tmp[tmp < 0] <- 0 + tmp + } else NA + # return return(list(aligned_image = aligned_image, + alignment_metrics = results[[1]], + matte_map = matte_map, + overlay_image = overlay_image, transformation = list( tfx_points = tfx_points, tfx_image = tfx_image @@ -3747,7 +4039,8 @@ getNonInteractiveRegistration <- function( channel_names = channel_names, query_ind = i, ref_ind = centre, - input = mapping_parameters + input = mapping_parameters, + compute_matte_map = FALSE ) } else { flag <- checkKeypoints(mapping_parameters$keypoints) @@ -3756,7 +4049,8 @@ getNonInteractiveRegistration <- function( keypoints_list = mapping_parameters$keypoints, query_ind = i, ref_ind = centre, - input = mapping_parameters + input = mapping_parameters, + compute_matte_map = FALSE ) } @@ -3781,3 +4075,59 @@ getNonInteractiveRegistration <- function( ) ) } + +#### +# Accuracy #### +#### + +#' getAlignmentAccuracy +#' +#' get accuracy measurements from two aligned images +#' +#' @param ref_image reference image +#' @param query_image query image +#' @param mask alignment mask +#' +#' @importFrom DelayedArray realize +#' @importFrom magick image_data +#' +#' @noRd +getAlignmentAccuracy <- function(ref_image, + query_image, + mask, + type, + compute_matte_map = TRUE){ + + # image info + ref_info <- getImageInfo(ref_image) + + # ref image + if (inherits(ref_image, "ImageArray")) { + ref_image <- DelayedArray::realize(ref_image) + ref_image <- array(as.raw(ref_image), dim = dim(ref_image)) + } else { + ref_image <- magick::image_data(ref_image, channels = "rgb") + } + + # query image + if (inherits(query_image, "ImageArray")) { + query_image <- DelayedArray::realize(query_image) + query_image <- array(as.raw(query_image), dim = dim(query_image)) + } else { + query_image <- magick::image_data(query_image, channels = "rgb") + } + + # mask + mask <- magick::image_data(mask, channels = "rgb") + + # calculate alignment accuracy + accuracy_rawvector(ref_image, + query_image, + mask, + width = ref_info$width, + height = ref_info$height, + type, + overlay_images = TRUE, + compute_matte_map = compute_matte_map) +} + diff --git a/docs/voltronobjects.html b/docs/voltronobjects.html index e9047c73..4b1e4b0b 100644 --- a/docs/voltronobjects.html +++ b/docs/voltronobjects.html @@ -683,8 +683,8 @@

Channel Names

vrImageChannelNames(melc_data)
 vrImageChannelNames(melc_data, assay = "MELC")
 vrImageChannelNames(melc_data, assay = "Assay1")
-
##        Assay    Layer         Sample Spatial  Channels
-## Assay1  MELC Section1 control_case_3    MELC DAPI,CD45
+
##        Assay    Layer         Sample Spatial           Channels
+## Assay1  MELC Section1 control_case_3    MELC DAPI,CD45,combined


@@ -726,6 +726,8 @@

Combining Image Channels

melc_data <- combineChannels(melc_data, 
                              channels = c("DAPI", "CD45"), colors = c("grey", "green"), 
                              channel_key = "combined")
+
## Warning in .local(object, ..., value = value): A channel with name 'combined' already exists in this vrImage object. 
+##  Overwriting ...

These new images can be stored as new channels within the same image object, and called later again

vrImageChannelNames(melc_data)
@@ -833,14 +835,14 @@

features

href="#assays">vrMainAssay)

selected_features <- vrFeatures(visium_data)
 selected_features[1:20]
-
##  [1] "Xkr4"          "Gm1992"        "Gm19938"       "Gm37381"       "Rp1"           "Sox17"         "Gm37587"       "Gm37323"       "Mrpl15"       
-## [10] "Lypla1"        "Tcea1"         "Rgs20"         "Gm16041"       "Atp6v1h"       "Oprk1"         "Npbwr1"        "Rb1cc1"        "4732440D04Rik"
-## [19] "Alkal1"        "St18"
+
##  [1] "Xkr4"          "Gm1992"        "Gm19938"       "Gm37381"       "Rp1"           "Sox17"         "Gm37587"      
+##  [8] "Gm37323"       "Mrpl15"        "Lypla1"        "Tcea1"         "Rgs20"         "Gm16041"       "Atp6v1h"      
+## [15] "Oprk1"         "Npbwr1"        "Rb1cc1"        "4732440D04Rik" "Alkal1"        "St18"
visium_data_subset <- subset(visium_data, features = selected_features[1:20])
 vrFeatures(visium_data_subset)
-
##  [1] "Xkr4"          "Gm1992"        "Gm19938"       "Gm37381"       "Rp1"           "Sox17"         "Gm37587"       "Gm37323"       "Mrpl15"       
-## [10] "Lypla1"        "Tcea1"         "Rgs20"         "Gm16041"       "Atp6v1h"       "Oprk1"         "Npbwr1"        "Rb1cc1"        "4732440D04Rik"
-## [19] "Alkal1"        "St18"
+
##  [1] "Xkr4"          "Gm1992"        "Gm19938"       "Gm37381"       "Rp1"           "Sox17"         "Gm37587"      
+##  [8] "Gm37323"       "Mrpl15"        "Lypla1"        "Tcea1"         "Rgs20"         "Gm16041"       "Atp6v1h"      
+## [15] "Oprk1"         "Npbwr1"        "Rb1cc1"        "4732440D04Rik" "Alkal1"        "St18"


diff --git a/inst/extdata/bspline_map.txt b/inst/extdata/bspline_map.txt index 12340a04..4b706fa8 100644 --- a/inst/extdata/bspline_map.txt +++ b/inst/extdata/bspline_map.txt @@ -34,4 +34,5 @@ (WriteResultImage "true") (WriteTransformParametersEachResolution "true") (ResultImageFormat "tif") -(ResultImagePixelType "unsigned char") \ No newline at end of file +(ResultImagePixelType "unsigned char") +(RandomSeed 121212) \ No newline at end of file diff --git a/man/getRcppManualRegistration.Rd b/man/getRcppManualRegistration.Rd index 9cb3fc8a..763e601f 100644 --- a/man/getRcppManualRegistration.Rd +++ b/man/getRcppManualRegistration.Rd @@ -9,6 +9,8 @@ getRcppManualRegistration( ref_image, query_landmark, reference_landmark, + invert_query = FALSE, + invert_ref = FALSE, method = "Homography", nonrigid = "TPS (OpenCV)" ) diff --git a/man/warpSimpleITKImage.Rd b/man/warpSimpleITKImage.Rd index 03bbddd9..68fd3d83 100644 --- a/man/warpSimpleITKImage.Rd +++ b/man/warpSimpleITKImage.Rd @@ -4,11 +4,9 @@ \alias{warpSimpleITKImage} \title{getRcppWarpImage} \usage{ -warpSimpleITKImage(ref_image, query_image, mapping) +warpSimpleITKImage(query_image, mapping) } \arguments{ -\item{ref_image}{reference image} - \item{query_image}{query image} \item{mapping}{a list of the homography matrices and TPS keypoints} diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 387e6cdb..20f8e9c6 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -11,9 +11,27 @@ Rcpp::Rostream& Rcpp::Rcout = Rcpp::Rcpp_cout_get(); Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); #endif +// accuracy_rawvector +Rcpp::List accuracy_rawvector(Rcpp::RawVector& ref_image, Rcpp::RawVector& query_image, Rcpp::RawVector& mask, const int width, const int height, std::string type, bool overlay_images, const bool compute_matte_map); +RcppExport SEXP _VoltRon_accuracy_rawvector(SEXP ref_imageSEXP, SEXP query_imageSEXP, SEXP maskSEXP, SEXP widthSEXP, SEXP heightSEXP, SEXP typeSEXP, SEXP overlay_imagesSEXP, SEXP compute_matte_mapSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< Rcpp::RawVector& >::type ref_image(ref_imageSEXP); + Rcpp::traits::input_parameter< Rcpp::RawVector& >::type query_image(query_imageSEXP); + Rcpp::traits::input_parameter< Rcpp::RawVector& >::type mask(maskSEXP); + Rcpp::traits::input_parameter< const int >::type width(widthSEXP); + Rcpp::traits::input_parameter< const int >::type height(heightSEXP); + Rcpp::traits::input_parameter< std::string >::type type(typeSEXP); + Rcpp::traits::input_parameter< bool >::type overlay_images(overlay_imagesSEXP); + Rcpp::traits::input_parameter< const bool >::type compute_matte_map(compute_matte_mapSEXP); + rcpp_result_gen = Rcpp::wrap(accuracy_rawvector(ref_image, query_image, mask, width, height, type, overlay_images, compute_matte_map)); + return rcpp_result_gen; +END_RCPP +} // automated_registeration_rawvector -Rcpp::List automated_registeration_rawvector(Rcpp::RawVector& ref_image, Rcpp::RawVector& query_image, const int width1, const int height1, const int width2, const int height2, const float GOOD_MATCH_PERCENT, const int MAX_FEATURES, const bool invert_query, const bool invert_ref, Rcpp::String flipflop_query, Rcpp::String flipflop_ref, Rcpp::String rotate_query, Rcpp::String rotate_ref, Rcpp::String matcher, Rcpp::String method, Rcpp::String nonrigid); -RcppExport SEXP _VoltRon_automated_registeration_rawvector(SEXP ref_imageSEXP, SEXP query_imageSEXP, SEXP width1SEXP, SEXP height1SEXP, SEXP width2SEXP, SEXP height2SEXP, SEXP GOOD_MATCH_PERCENTSEXP, SEXP MAX_FEATURESSEXP, SEXP invert_querySEXP, SEXP invert_refSEXP, SEXP flipflop_querySEXP, SEXP flipflop_refSEXP, SEXP rotate_querySEXP, SEXP rotate_refSEXP, SEXP matcherSEXP, SEXP methodSEXP, SEXP nonrigidSEXP) { +Rcpp::List automated_registeration_rawvector(Rcpp::RawVector& ref_image, Rcpp::RawVector& query_image, const int width1, const int height1, const int width2, const int height2, const float GOOD_MATCH_PERCENT, const int MAX_FEATURES, const bool invert_query, const bool invert_ref, Rcpp::String flipflop_query, Rcpp::String flipflop_ref, Rcpp::String rotate_query, Rcpp::String rotate_ref, Rcpp::String matcher, Rcpp::String method, Rcpp::String nonrigid, const bool compute_matte_map); +RcppExport SEXP _VoltRon_automated_registeration_rawvector(SEXP ref_imageSEXP, SEXP query_imageSEXP, SEXP width1SEXP, SEXP height1SEXP, SEXP width2SEXP, SEXP height2SEXP, SEXP GOOD_MATCH_PERCENTSEXP, SEXP MAX_FEATURESSEXP, SEXP invert_querySEXP, SEXP invert_refSEXP, SEXP flipflop_querySEXP, SEXP flipflop_refSEXP, SEXP rotate_querySEXP, SEXP rotate_refSEXP, SEXP matcherSEXP, SEXP methodSEXP, SEXP nonrigidSEXP, SEXP compute_matte_mapSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -34,7 +52,8 @@ BEGIN_RCPP Rcpp::traits::input_parameter< Rcpp::String >::type matcher(matcherSEXP); Rcpp::traits::input_parameter< Rcpp::String >::type method(methodSEXP); Rcpp::traits::input_parameter< Rcpp::String >::type nonrigid(nonrigidSEXP); - rcpp_result_gen = Rcpp::wrap(automated_registeration_rawvector(ref_image, query_image, width1, height1, width2, height2, GOOD_MATCH_PERCENT, MAX_FEATURES, invert_query, invert_ref, flipflop_query, flipflop_ref, rotate_query, rotate_ref, matcher, method, nonrigid)); + Rcpp::traits::input_parameter< const bool >::type compute_matte_map(compute_matte_mapSEXP); + rcpp_result_gen = Rcpp::wrap(automated_registeration_rawvector(ref_image, query_image, width1, height1, width2, height2, GOOD_MATCH_PERCENT, MAX_FEATURES, invert_query, invert_ref, flipflop_query, flipflop_ref, rotate_query, rotate_ref, matcher, method, nonrigid, compute_matte_map)); return rcpp_result_gen; END_RCPP } @@ -101,8 +120,8 @@ BEGIN_RCPP END_RCPP } // manual_registeration_rawvector -Rcpp::List manual_registeration_rawvector(Rcpp::RawVector ref_image, Rcpp::RawVector query_image, Rcpp::NumericMatrix reference_landmark, Rcpp::NumericMatrix query_landmark, const int width1, const int height1, const int width2, const int height2, Rcpp::String method, Rcpp::String nonrigid); -RcppExport SEXP _VoltRon_manual_registeration_rawvector(SEXP ref_imageSEXP, SEXP query_imageSEXP, SEXP reference_landmarkSEXP, SEXP query_landmarkSEXP, SEXP width1SEXP, SEXP height1SEXP, SEXP width2SEXP, SEXP height2SEXP, SEXP methodSEXP, SEXP nonrigidSEXP) { +Rcpp::List manual_registeration_rawvector(Rcpp::RawVector ref_image, Rcpp::RawVector query_image, Rcpp::NumericMatrix reference_landmark, Rcpp::NumericMatrix query_landmark, const int width1, const int height1, const int width2, const int height2, const bool invert_query, const bool invert_ref, Rcpp::String method, Rcpp::String nonrigid, const bool compute_matte_map); +RcppExport SEXP _VoltRon_manual_registeration_rawvector(SEXP ref_imageSEXP, SEXP query_imageSEXP, SEXP reference_landmarkSEXP, SEXP query_landmarkSEXP, SEXP width1SEXP, SEXP height1SEXP, SEXP width2SEXP, SEXP height2SEXP, SEXP invert_querySEXP, SEXP invert_refSEXP, SEXP methodSEXP, SEXP nonrigidSEXP, SEXP compute_matte_mapSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -114,9 +133,12 @@ BEGIN_RCPP Rcpp::traits::input_parameter< const int >::type height1(height1SEXP); Rcpp::traits::input_parameter< const int >::type width2(width2SEXP); Rcpp::traits::input_parameter< const int >::type height2(height2SEXP); + Rcpp::traits::input_parameter< const bool >::type invert_query(invert_querySEXP); + Rcpp::traits::input_parameter< const bool >::type invert_ref(invert_refSEXP); Rcpp::traits::input_parameter< Rcpp::String >::type method(methodSEXP); Rcpp::traits::input_parameter< Rcpp::String >::type nonrigid(nonrigidSEXP); - rcpp_result_gen = Rcpp::wrap(manual_registeration_rawvector(ref_image, query_image, reference_landmark, query_landmark, width1, height1, width2, height2, method, nonrigid)); + Rcpp::traits::input_parameter< const bool >::type compute_matte_map(compute_matte_mapSEXP); + rcpp_result_gen = Rcpp::wrap(manual_registeration_rawvector(ref_image, query_image, reference_landmark, query_landmark, width1, height1, width2, height2, invert_query, invert_ref, method, nonrigid, compute_matte_map)); return rcpp_result_gen; END_RCPP } @@ -147,6 +169,19 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// generateOverlapMask +Rcpp::IntegerVector generateOverlapMask(Rcpp::NumericVector& dsize, Rcpp::NumericMatrix& trans_mat, Rcpp::NumericVector& ssize); +RcppExport SEXP _VoltRon_generateOverlapMask(SEXP dsizeSEXP, SEXP trans_matSEXP, SEXP ssizeSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< Rcpp::NumericVector& >::type dsize(dsizeSEXP); + Rcpp::traits::input_parameter< Rcpp::NumericMatrix& >::type trans_mat(trans_matSEXP); + Rcpp::traits::input_parameter< Rcpp::NumericVector& >::type ssize(ssizeSEXP); + rcpp_result_gen = Rcpp::wrap(generateOverlapMask(dsize, trans_mat, ssize)); + return rcpp_result_gen; +END_RCPP +} // build_snn_rank Rcpp::List build_snn_rank(Rcpp::IntegerMatrix neighbors); RcppExport SEXP _VoltRon_build_snn_rank(SEXP neighborsSEXP) { @@ -182,14 +217,16 @@ END_RCPP } static const R_CallMethodDef CallEntries[] = { - {"_VoltRon_automated_registeration_rawvector", (DL_FUNC) &_VoltRon_automated_registeration_rawvector, 17}, + {"_VoltRon_accuracy_rawvector", (DL_FUNC) &_VoltRon_accuracy_rawvector, 8}, + {"_VoltRon_automated_registeration_rawvector", (DL_FUNC) &_VoltRon_automated_registeration_rawvector, 18}, {"_VoltRon_replaceNaMatrix", (DL_FUNC) &_VoltRon_replaceNaMatrix, 2}, {"_VoltRon_warpRcppImage", (DL_FUNC) &_VoltRon_warpRcppImage, 7}, {"_VoltRon_warpImageAuto", (DL_FUNC) &_VoltRon_warpImageAuto, 7}, {"_VoltRon_warpImageManual", (DL_FUNC) &_VoltRon_warpImageManual, 7}, - {"_VoltRon_manual_registeration_rawvector", (DL_FUNC) &_VoltRon_manual_registeration_rawvector, 10}, + {"_VoltRon_manual_registeration_rawvector", (DL_FUNC) &_VoltRon_manual_registeration_rawvector, 13}, {"_VoltRon_manual_registeration_matrix", (DL_FUNC) &_VoltRon_manual_registeration_matrix, 5}, {"_VoltRon_applyRcppMapping", (DL_FUNC) &_VoltRon_applyRcppMapping, 2}, + {"_VoltRon_generateOverlapMask", (DL_FUNC) &_VoltRon_generateOverlapMask, 3}, {"_VoltRon_build_snn_rank", (DL_FUNC) &_VoltRon_build_snn_rank, 1}, {"_VoltRon_build_snn_number", (DL_FUNC) &_VoltRon_build_snn_number, 1}, {"_VoltRon_replacePatternInRcppVectorWrapper", (DL_FUNC) &_VoltRon_replacePatternInRcppVectorWrapper, 3}, diff --git a/src/accuracy.cpp b/src/accuracy.cpp new file mode 100644 index 00000000..1374d4fa --- /dev/null +++ b/src/accuracy.cpp @@ -0,0 +1,66 @@ +#include + +// OpenCV +#include + +// Library +#include "auxiliary.h" +#include "image.h" +#include "metrics.h" +#include "matte_mi.h" + +// Namespaces +using namespace Rcpp; +using namespace std; +using namespace cv; + +// [[Rcpp::export]] +Rcpp::List accuracy_rawvector(Rcpp::RawVector& ref_image, + Rcpp::RawVector& query_image, + Rcpp::RawVector& mask, + const int width, + const int height, + std::string type, + bool overlay_images = true, + const bool compute_matte_map = true) { + // results + Rcpp::List out(3); + + // Read images + cv::Mat imReference = imageToMat(ref_image, width, height); + cv::Mat imReg = imageToMat(query_image, width, height); + cv::Mat maskReg = imageToMat(mask, width, height); + + // process + Mat im1Proc, im2Proc; + cvtColor(imReg, im1Proc, cv::COLOR_BGR2GRAY); + cvtColor(imReference, im2Proc, cv::COLOR_BGR2GRAY); + cvtColor(maskReg, maskReg, cv::COLOR_BGR2GRAY); + + // get metrics + std::map accuracy; + accuracy = getAlignmentMetrics(im1Proc, im2Proc, maskReg, type); + out[0] = accuracy; + + // get matte map + Mat1d accuracyMatte; + if(compute_matte_map){ + accuracyMatte = MatteMIMap(im2Proc, im1Proc, maskReg, 50); + out[1] = matToNumericMatrix(accuracyMatte); // Matte MI metric + } else { + out[1] = R_NilValue; + } + + // image overlay + if(overlay_images){ + cv::addWeighted(im2Proc, 0.7, im1Proc, 0.3, 0, im1Proc); + cvtColor(im1Proc, im1Proc, cv::COLOR_GRAY2BGR); + im1Proc = resize_image(im1Proc, 500); + out[2] = matToImage(im1Proc); + } else { + out[2] = R_NilValue; + } + + // return + return out; +} \ No newline at end of file diff --git a/src/automated_registration.cpp b/src/automated_registration.cpp index 7bbb1305..dfec93c0 100644 --- a/src/automated_registration.cpp +++ b/src/automated_registration.cpp @@ -4,17 +4,22 @@ #include #include "opencv2/features2d.hpp" #include "opencv2/shape/shape_transformer.hpp" -// #include -// Internal functions +// Library #include "auxiliary.h" #include "image.h" +#include "metrics.h" +#include "matte_mi.h" // Namespaces using namespace Rcpp; using namespace std; using namespace cv; +//// +// Parameters +//// + // SIFT Parameters struct SIFTParameters { @@ -27,95 +32,9 @@ struct SIFTParameters const int ransac_maxIters=2000; }; -// check if keypoints are degenerate -bool check_degenerate(std::vector &points1, std::vector &points2) { - - // get sd - double points1_sd = cppSD(points1); - double points2_sd = cppSD(points2); - - // get warning message - bool is_degenerate = FALSE; - if(points1_sd < 1.0 | points2_sd < 1.0){ - is_degenerate = TRUE; - Rcout << "WARNING: points may be in a degenerate configuration." << endl; - } - - return is_degenerate; -} - -// check distribution of registered points -std::string check_transformation_by_point_distribution(Mat &im, Mat &h){ - - // message - std::string message; - - // get image shape - int height = im.rows; - int width = im.cols; - int height_interval = height > 50 ? (double) height/50.0 : 1; - int width_interval = width > 50 ? (double) width/50.0 : 1; - - // perspective transformation of grid points - std::vector gridpoints; - for (double i = 0.0; i <= height; i += height_interval) { - for (double j = 0.0; j <= width; j += width_interval) { - gridpoints.push_back(cv::Point2f(j,i)); - } - } - - // register grid points - std::vector gridpoints_reg; - if (h.rows == 2){ - cv::transform(gridpoints, gridpoints_reg, h); - } else if(h.rows == 3) { - cv::perspectiveTransform(gridpoints, gridpoints_reg, h); - } else { - message = "no distribution"; - return message; - } - - // Compute the standard deviation of the transformed points - double gridpoints_reg_sd = cppSD(gridpoints_reg); - - // get warning message - if(gridpoints_reg_sd < 1.0 | gridpoints_reg_sd > max(height, width)){ - message = "large distribution"; - Rcout << "WARNING: Transformation may be poor - transformed points grid seem to be concentrated!" << endl; - } else { - message = "small distribution"; - } - - return message; -} - -bool check_matches(Mat &mask){ - int j=0; - for (int i = 0; i < mask.rows; i++) { - if (mask.at(i)) { - j++; - } - } - return j > 6; -} - -// do overall checks on keypoints and images -bool check_transformation_metrics(std::vector &points1, std::vector &points2, Mat &im2, Mat &h, Mat &mask) { - - // check keypoint standard deviation - bool is_degenerate = check_degenerate(points1, points2); - - // TODO: check transformation - // make keypoints from points - // std::string transformation; - // transformation = check_transformation_by_pts_mean_sqrt(keypoints1, keypoints2, h, mask); - - // check distribution - std::string distribution; - distribution = check_transformation_by_point_distribution(im2, h); - - return is_degenerate; -} +//// +// Manage Keypoints and Matches +//// // get good matching keypoints void getGoodMatches(std::vector> &matches12,std::vector> &matches21, @@ -171,17 +90,10 @@ void getGoodMatches(std::vector> &matches12,std::vector, std::vector>().swap(matches21_map); } -void getGoodMatches_temp(std::vector> matches, std::vector &good_matches, const float lowe_ratio = 0.8) -{ - for (size_t i = 0; i < matches.size(); i++) { - if (matches[i][0].distance < lowe_ratio * matches[i][1].distance) { - good_matches.push_back(matches[i][0]); - } - } -} - // remove duplicate keypoints for TPS -void removeCloseMatches(std::vector& points1, std::vector& points2, float threshold = std::numeric_limits::epsilon()) { +void removeCloseMatches(std::vector& points1, + std::vector& points2, + float threshold = std::numeric_limits::epsilon()) { // Create a vector to store filtered points std::vector filtered_points1; @@ -293,6 +205,10 @@ void keepTopKeypoints(std::vector &keypoints, Mat &descriptors, SIFTPa } } +//// +// Compute SIFT/ORB and transformations +//// + void computeSIFTTiles(Mat &im, std::vector &keypoints, Mat &descriptors, Ptr &sift, SIFTParameters params){ // profiler @@ -352,6 +268,7 @@ void computeSIFTTiles(Mat &im, std::vector &keypoints, Mat &descriptor } } + bool getSIFTTransformationMatrixSingle( Mat &im1Proc, Mat &im2Proc, Mat &h, Mat &mask, Mat &imMatches, @@ -370,8 +287,6 @@ bool getSIFTTransformationMatrixSingle( Ptr sift = cv::SIFT::create(params.sift_nfeatures); computeSIFTTiles(im1Proc, keypoints1, descriptors1, sift, params); computeSIFTTiles(im2Proc, keypoints2, descriptors2, sift, params); - Rcout << "MESSAGE: Generated " << keypoints1.size() << " and " << keypoints2.size() << " keypoints" << endl; - Rcout << "DONE: SIFT based key-points detection and descriptors computation" << endl; // filter duplicates filterDuplicateKeypoints(keypoints1, descriptors1); @@ -380,7 +295,6 @@ bool getSIFTTransformationMatrixSingle( // get top key points keepTopKeypoints(keypoints1, descriptors1, params); keepTopKeypoints(keypoints2, descriptors2, params); - Rcout << "MESSAGE: Filtered other than " << keypoints1.size() << " and " << keypoints2.size() << " keypoints" << endl; /////////////////////// /// Compute FLANN ///// @@ -389,7 +303,6 @@ bool getSIFTTransformationMatrixSingle( // Match features using FLANN matching std::vector> matches12, matches21; getFLANNMatches(descriptors1, descriptors2, matches12, matches21); - Rcout << "DONE: FLANN - Fast Library for Approximate Nearest Neighbors - descriptor matching" << endl; // TODO: can I release there now ? descriptors1.release(); @@ -398,12 +311,7 @@ bool getSIFTTransformationMatrixSingle( // Find good matches std::vector good_matches; getGoodMatches(matches12, matches21, good_matches); - Rcout << "DONE: get good mutual matches by distance thresholding" << endl; - // TODO: can I release there now ? - std::vector>().swap(matches12); - std::vector>().swap(matches21); - /////////////////////// /// Find Homography /// /////////////////////// @@ -414,10 +322,7 @@ bool getSIFTTransformationMatrixSingle( points1.push_back(keypoints1[good_matches[i].queryIdx].pt); points2.push_back(keypoints2[good_matches[i].trainIdx].pt); } - - // check variable - Rcout << "MESSAGE: Calculating" << (run_Affine ? " (Affine) " : " (Homography) ") << "Transformation Matrix" << endl; - + // Find transformation matrix if(points1.size() > 0){ if(run_Affine){ @@ -444,35 +349,31 @@ bool getSIFTTransformationMatrixSingle( return false; } - // Draw top matches and good ones only - std::vector top_matches; - std::vector keypoints1_best, keypoints2_best; + // filter keypoints by good matches + std::vector keypoints1_good, keypoints2_good; for(size_t i = 0; i < good_matches.size(); i++ ) { - keypoints1_best.push_back(keypoints1[good_matches[i].queryIdx]); - keypoints2_best.push_back(keypoints2[good_matches[i].trainIdx]); - } - std::vector keypoints1_best2, keypoints2_best2; - int j=0; - for (int i = 0; i < mask.rows; i++) { - if (mask.at(i)) { - keypoints1_best2.push_back(keypoints1_best[i]); - keypoints2_best2.push_back(keypoints2_best[i]); - top_matches.push_back(cv::DMatch(static_cast(j), static_cast(j), 0)); - j++; - } + keypoints1_good.push_back(keypoints1[good_matches[i].queryIdx]); + keypoints2_good.push_back(keypoints2[good_matches[i].trainIdx]); } - scaledDrawMatches(im1Proc, keypoints1_best2, im2Proc, keypoints2_best2, top_matches, imMatches); - - // TODO: can I release there now ? - // std::vector().swap(keypoints1_best); - // std::vector().swap(keypoints2_best); - // std::vector().swap(keypoints1_best2); - // std::vector().swap(keypoints2_best2); - // std::vector().swap(top_matches); + + // filter keypoints using mask + std::vector keypoints1_masked, keypoints2_masked; + std::vector top_matches; + maskKeypoints(keypoints1_good, keypoints2_good, + keypoints1_masked, keypoints2_masked, + top_matches, mask); + + // convert keypoints to points + points1 = KeyPointToPoint2f(keypoints1_masked); + points2 = KeyPointToPoint2f(keypoints2_masked); + + // draw matches + scaledDrawMatches(im1Proc, keypoints1_masked, im2Proc, keypoints2_masked, + top_matches, imMatches); // check number of matches - return check_matches(mask); + return checkMaskAbundance(mask); } void getSIFTTransformationMatrix( @@ -493,67 +394,73 @@ void getSIFTTransformationMatrix( // check variable bool check; - Rcout << "MESSAGE: Calculating" << (run_Affine ? " (Affine) " : " (Homography) ") << "Transformation Matrix" << endl; + Rcout << "Calculating" << (run_Affine ? " (Affine) " : " (Homography) ") << "Transformation Matrix" << endl; + + Rcout << "Round 1: No histogram equalization" << endl; // find matches and points check = getSIFTTransformationMatrixSingle(im1Proc, im2Proc, h, mask, imMatches, points1, points2, run_Affine, params, is_faulty); - Rcout << "DONE: calculated homography matrix with " << points1.size() << " points" << endl; - + // equalize first image if fails if(!check){ + + // clear points, mask is reset itself + points1.clear(); + points2.clear(); Mat im1Proc_eq; cv::equalizeHist(im1Proc, im1Proc_eq); - Rcout << "MESSAGE: Calculating Transformation Matrix with histogram equalization (1)" << endl; + Rcout << "Round 2: Histogram equalization of Image 1" << endl; check = getSIFTTransformationMatrixSingle(im1Proc_eq, im2Proc, h, mask, imMatches, points1, points2, run_Affine, params, is_faulty); - Rcout << "DONE: calculated homography matrix with " << points1.size() << " points" << endl; } else { return; } // equalize second image if fails if(!check){ + + // clear points, mask is reset itself + points1.clear(); + points2.clear(); cv::equalizeHist(im2Proc, im2Proc_eq); - Rcout << "MESSAGE: Calculating Transformation Matrix with histogram equalization (2)" << endl; + Rcout << "Round 3: Histogram equalization of Image 2" << endl; check = getSIFTTransformationMatrixSingle(im1Proc, im2Proc_eq, h, mask, imMatches, points1, points2, run_Affine, params, is_faulty); - Rcout << "DONE: calculated homography matrix with " << points1.size() << " points" << endl; } else { return; } // last try with both equalized images if(!check){ + + // clear points, mask is reset itself + points1.clear(); + points2.clear(); cv::equalizeHist(im1Proc, im1Proc_eq2); cv::equalizeHist(im2Proc, im2Proc_eq2); - Rcout << "MESSAGE: Calculating Transformation Matrix with histogram equalization (3)" << endl; + Rcout << "Round 4: Histogram equalization of Image 1 and 2" << endl; check = getSIFTTransformationMatrixSingle(im1Proc_eq2, im2Proc_eq2, h, mask, imMatches, points1, points2, run_Affine, params, is_faulty); - Rcout << "DONE: calculated homography matrix with " << points1.size() << " points" << endl; } else { return; } - // TODO: release ? - // im1Proc_eq.release(); - // im1Proc_eq2.release(); - // im2Proc_eq.release(); - // im2Proc_eq2.release(); + } bool getORBTransformationMatrix( @@ -573,22 +480,19 @@ bool getORBTransformationMatrix( Ptr orb = ORB::create(MAX_FEATURES); orb->detectAndCompute(im1Proc, Mat(), keypoints1, descriptors1); orb->detectAndCompute(im2Proc, Mat(), keypoints2, descriptors2); - Rcout << "DONE: orb based key-points detection and descriptors computation" << endl; - + // Match features. std::vector matches; Ptr matcher = DescriptorMatcher::create("BruteForce-Hamming"); matcher->match(descriptors1, descriptors2, matches, Mat()); - Rcout << "DONE: BruteForce-Hamming - descriptor matching" << endl; - + // Sort matches by score std::sort(matches.begin(), matches.end()); // Remove not so good matches const int numGoodMatches = matches.size() * GOOD_MATCH_PERCENT; matches.erase(matches.begin()+numGoodMatches, matches.end()); - Rcout << "DONE: get good matches by distance thresholding" << endl; - + // Extract location of good matches for( size_t i = 0; i < matches.size(); i++ ) { @@ -597,8 +501,8 @@ bool getORBTransformationMatrix( } // check variable - Rcout << "MESSAGE: Calculating" << (run_Affine ? " (Affine) " : " (Homography) ") << "Transformation Matrix" << endl; - + Rcout << "Calculating" << (run_Affine ? " (Affine) " : " (Homography) ") << "Transformation Matrix" << endl; + // Find transformation matrix if(points1.size() > 0){ if(run_Affine){ @@ -616,43 +520,64 @@ bool getORBTransformationMatrix( mask); } } else { - Rcout << "Found no matches!" << endl; + Rcout << "WARNING: Found no matches!" << endl; return false; } - - // Draw top matches and good ones only - std::vector top_matches; - std::vector keypoints1_best, keypoints2_best; + + // filter keypoints by good matches + std::vector keypoints1_good, keypoints2_good; for(size_t i = 0; i < matches.size(); i++ ) { - keypoints1_best.push_back(keypoints1[matches[i].queryIdx]); - keypoints2_best.push_back(keypoints2[matches[i].trainIdx]); + keypoints1_good.push_back(keypoints1[matches[i].queryIdx]); + keypoints2_good.push_back(keypoints2[matches[i].trainIdx]); } - std::vector keypoints1_best2, keypoints2_best2; - int j=0; - for (int i = 0; i < mask.rows; i++) { - if (mask.at(i)) { - keypoints1_best2.push_back(keypoints1_best[i]); - keypoints2_best2.push_back(keypoints2_best[i]); - top_matches.push_back(cv::DMatch(static_cast(j), static_cast(j), 0)); - j++; - } - } - scaledDrawMatches(im1Proc, keypoints1_best2, im2Proc, keypoints2_best2, top_matches, imMatches); - + + // filter keypoints using mask + std::vector keypoints1_masked, keypoints2_masked; + std::vector top_matches; + maskKeypoints(keypoints1_good, keypoints2_good, + keypoints1_masked, keypoints2_masked, + top_matches, mask); + + // convert keypoints to points + points1 = KeyPointToPoint2f(keypoints1_masked); + points2 = KeyPointToPoint2f(keypoints2_masked); + + // draw matches + scaledDrawMatches(im1Proc, keypoints1_masked, im2Proc, keypoints2_masked, + top_matches, imMatches); + // check number of matches - return check_matches(mask); + return checkMaskAbundance(mask); } +//// +// Align Images +//// + // align images with FLANN algorithm -void alignImages(Mat &im1, Mat &im2, Mat &im1Reg, Mat &im1Overlay, - Mat &imMatches, Mat &h, Rcpp::List &keypoints, - const float GOOD_MATCH_PERCENT, const int MAX_FEATURES, +void alignImages(Mat &im1, + Mat &im2, + Mat &im1Reg, + Mat &im1Overlay, + Mat &imMatches, + Mat &h, + Rcpp::List &keypoints, + const float GOOD_MATCH_PERCENT, + const int MAX_FEATURES, Rcpp::String matcher, - const bool invert_query, const bool invert_ref, - const char* flipflop_query, const char* flipflop_ref, - const char* rotate_query, const char* rotate_ref, - const bool run_Affine, const bool run_TPS) + const bool invert_query, + const bool invert_ref, + const char* flipflop_query, + const char* flipflop_ref, + const char* rotate_query, + const char* rotate_ref, + const bool run_Affine, + const bool run_TPS, + const bool compute_matte_map, + Mat1d &accuracyMatte, + std::map &accuracy_coarse, + std::map &accuracy_fine) { // parameters @@ -683,7 +608,7 @@ void alignImages(Mat &im1, Mat &im2, Mat &im1Reg, Mat &im1Overlay, if(strcmp(matcher.get_cstring(), "BRUTE-FORCE") == 0){ // message - Rcout << "MESSAGE: Running BRUTE-FORCE Alignment" << endl; + Rcout << "Running Coarse Alignment (BRUTE-FORCE)" << endl; // run ORB bool check; @@ -694,32 +619,57 @@ void alignImages(Mat &im1, Mat &im2, Mat &im1Reg, Mat &im1Overlay, } else { // message - Rcout << "MESSAGE: Running SIFT+FLANN Alignment" << ((run_TPS) ? " with TPS" : "") << endl; + // Rcout << "Running SIFT+FLANN Alignment" << ((run_TPS) ? " with TPS" : "") << endl; + Rcout << "Running Coarse Alignment (SIFT+FLANN)" << endl; // run SIFT getSIFTTransformationMatrix(im1Proc, im2Proc, h, mask, imMatches, points1, points2, run_Affine, is_faulty); } - - // check result - is_faulty = check_transformation_metrics(points1, points2, im2, h, mask); - Rcout << "MESSAGE: Registration is " << (is_faulty ? "degenerate!" : "not degenerate!") << endl; - // Use homography to warp image - if(h.rows == 2){ - warpAffine(im1Proc, im1Proc, h, im2Proc.size()); - warpAffine(im1NormalProc, im1NormalProc, h, im2Proc.size()); - } else if(h.rows == 3){ - warpPerspective(im1Proc, im1Proc, h, im2Proc.size()); - warpPerspective(im1NormalProc, im1NormalProc, h, im2Proc.size()); + // warp mask and image + cv::Mat alignmentMask; + std::map keypoint_metrics; + if(!h.empty()){ + alignmentMask = generateOverlapMask(im2Proc.size(), + h, + im1Proc.size()); + + // get keypoint metrics before warping + keypoint_metrics = getKeypointMetrics(points1, points2, + im1Proc, im2Proc, h, mask); + is_faulty = (bool) keypoint_metrics["Degenerate"]; + + if(h.rows == 2){ + warpAffine(im1Proc, im1Proc, h, im2Proc.size()); + warpAffine(im1NormalProc, im1NormalProc, h, im2Proc.size()); + } else if(h.rows == 3){ + warpPerspective(im1Proc, im1Proc, h, im2Proc.size()); + warpPerspective(im1NormalProc, im1NormalProc, h, im2Proc.size()); + } } else { Rcout << "WARNING: No transformation was found" << endl; return; } - - Rcout << "DONE: warped query image" << endl; + // get alignment metrics + std::map image_metrics; + image_metrics = getAlignmentMetrics(im1Proc, im2Proc, + alignmentMask, "Coarse"); + + // combine metrics + std::vector> temp_map; + temp_map.reserve(keypoint_metrics.size() + image_metrics.size()); + std::copy(keypoint_metrics.begin(), keypoint_metrics.end(), std::back_inserter(temp_map)); + std::copy(image_metrics.begin(), image_metrics.end(), std::back_inserter(temp_map)); + std::map final_map(temp_map.begin(), temp_map.end()); + accuracy_coarse = final_map; + + // get matte metric + if(compute_matte_map) + accuracyMatte = MatteMIMap(im2Proc, im1Proc, alignmentMask, 50); + /////////////////////// /// Find Homography /// /////////////////////// @@ -736,55 +686,47 @@ void alignImages(Mat &im1, Mat &im2, Mat &im1Reg, Mat &im1Overlay, // TPS is requested (only if FLANN succeeded) } else { - Rcout << "MESSAGE: Running Thin-Plate-Spline Alignment" << endl; - - // Filtered points (inliers) based on the mask - std::vector filtered_points1; - std::vector filtered_points2; - for (int i = 0; i < mask.rows; i++) { - if (mask.at(i)) { - filtered_points1.push_back(points1[i]); - filtered_points2.push_back(points2[i]); - } - } - removeCloseMatches(filtered_points1, filtered_points2); + Rcout << "Running Fine Alignment (Thin-Plate-Spline)" << endl; + // remove close looking matches + removeCloseMatches(points1, points2); + // transform query - std::vector filtered_points1_reg; + std::vector points1_reg; if (h.rows == 2){ - cv::transform(filtered_points1, filtered_points1_reg, h); + cv::transform(points1, points1_reg, h); } else { - cv::perspectiveTransform(filtered_points1, filtered_points1_reg, h); + cv::perspectiveTransform(points1, points1_reg, h); } // get TPS matches std::vector matches; - for (unsigned int i = 0; i < filtered_points2.size(); i++) + for (unsigned int i = 0; i < points2.size(); i++) matches.push_back(cv::DMatch(i, i, 0)); // calculate TPS transformation Ptr tps = cv::createThinPlateSplineShapeTransformer(0); - tps->estimateTransformation(filtered_points2, filtered_points1_reg, matches); + tps->estimateTransformation(points2, points1_reg, matches); // save keypoints - keypoints[0] = point2fToNumericMatrix(filtered_points2); - keypoints[1] = point2fToNumericMatrix(filtered_points1_reg); - - // determine extension limits for both images - int y_max = max(im1Proc.rows, im2.rows); - int x_max = max(im1Proc.cols, im2.cols); + keypoints[0] = point2fToNumericMatrix(points2); + keypoints[1] = point2fToNumericMatrix(points1_reg); - // extend images - cv::copyMakeBorder(im1Proc, im1Proc, 0.0, (int) (y_max - im1Proc.rows), 0.0, (x_max - im1Proc.cols), cv::BORDER_CONSTANT, Scalar(0, 0, 0)); - cv::copyMakeBorder(im1NormalProc, im1NormalProc, 0.0, (int) (y_max - im1NormalProc.rows), 0.0, (x_max - im1NormalProc.cols), cv::BORDER_CONSTANT, Scalar(0, 0, 0)); - - // transform image - tps->warpImage(im1Proc, im1Proc); - tps->warpImage(im1NormalProc, im1NormalProc); + // warp overlap mask + alignmentMask = warpTPSImage(im2Proc, alignmentMask, tps, + im2Proc.rows, im2Proc.cols, + cv::INTER_NEAREST); + + // transform image using trained tps + im1Proc = warpTPSImage(im2Proc, im1Proc, tps, + im2Proc.rows, im2Proc.cols, cv::INTER_LINEAR); + im1NormalProc = warpTPSImage(im2Proc, im1NormalProc, tps, + im2Proc.rows, im2Proc.cols, cv::INTER_LINEAR); - // resize image - im1Proc = im1Proc(cv::Range(0,im2Proc.size().height), cv::Range(0,im2Proc.size().width)); - im1NormalProc = im1NormalProc(cv::Range(0,im2Proc.size().height), cv::Range(0,im2Proc.size().width)); + // get matte metric, process + accuracy_fine = getAlignmentMetrics(im1Proc, im2Proc, alignmentMask, "Fine"); + if(compute_matte_map) + accuracyMatte = MatteMIMap(im2Proc, im1Proc, alignmentMask, 50); // change color map cv::addWeighted(im2Proc, 0.7, im1Proc, 0.3, 0, im1Proc); @@ -807,20 +749,32 @@ void alignImages(Mat &im1, Mat &im2, Mat &im1Reg, Mat &im1Overlay, } // [[Rcpp::export]] -Rcpp::List automated_registeration_rawvector(Rcpp::RawVector& ref_image, Rcpp::RawVector& query_image, - const int width1, const int height1, - const int width2, const int height2, - const float GOOD_MATCH_PERCENT, const int MAX_FEATURES, - const bool invert_query, const bool invert_ref, - Rcpp::String flipflop_query, Rcpp::String flipflop_ref, - Rcpp::String rotate_query, Rcpp::String rotate_ref, - Rcpp::String matcher, Rcpp::String method, Rcpp::String nonrigid) +Rcpp::List automated_registeration_rawvector(Rcpp::RawVector& ref_image, + Rcpp::RawVector& query_image, + const int width1, + const int height1, + const int width2, + const int height2, + const float GOOD_MATCH_PERCENT, + const int MAX_FEATURES, + const bool invert_query, + const bool invert_ref, + Rcpp::String flipflop_query, + Rcpp::String flipflop_ref, + Rcpp::String rotate_query, + Rcpp::String rotate_ref, + Rcpp::String matcher, + Rcpp::String method, + Rcpp::String nonrigid, + const bool compute_matte_map = true) { // Return data - Rcpp::List out(5); + Rcpp::List out(8); Rcpp::List out_trans(2); Rcpp::List keypoints(2); Mat imOverlay, imReg, h, imMatches; + Mat1d accuracyMatte; + std::map accuracy_coarse, accuracy_fine; // Read reference image cv::Mat imReference = imageToMat(ref_image, width1, height1); @@ -841,7 +795,11 @@ Rcpp::List automated_registeration_rawvector(Rcpp::RawVector& ref_image, Rcpp::R invert_query, invert_ref, flipflop_query.get_cstring(), flipflop_ref.get_cstring(), rotate_query.get_cstring(), rotate_ref.get_cstring(), - run_Affine, run_TPS); + run_Affine, run_TPS, + compute_matte_map, + accuracyMatte, + accuracy_coarse, + accuracy_fine); // transformation matrix, can be either a matrix, set of keypoints or both out_trans[0] = matToNumericMatrix(h.clone()); @@ -857,10 +815,16 @@ Rcpp::List automated_registeration_rawvector(Rcpp::RawVector& ref_image, Rcpp::R out[2] = matToImage(imReg); // registered image out[3] = matToImage(imMatches); // keypoint matching image out[4] = matToImage(imOverlay); // overlay image + out[5] = matToNumericMatrix(accuracyMatte); // Matte MI metric + out[6] = accuracy_coarse; // accuracy scores (coarse) + out[7] = accuracy_fine; // accuracy scores (fine) } else { out[2] = R_NilValue; out[3] = R_NilValue; out[4] = R_NilValue; + out[5] = R_NilValue; + out[6] = R_NilValue; + out[7] = R_NilValue; } // release @@ -872,498 +836,4 @@ Rcpp::List automated_registeration_rawvector(Rcpp::RawVector& ref_image, Rcpp::R // return return out; -} - -///////////////// -/// scratch ///// -///////////////// - -// align images with BRUTE FORCE algorithm -void alignImagesBRUTE(Mat &im1, Mat &im2, Mat &im1Reg, Mat &im1Overlay, Mat &imMatches, Mat &h, - const float GOOD_MATCH_PERCENT, const int MAX_FEATURES, - const bool invert_query, const bool invert_ref, - const char* flipflop_query, const char* flipflop_ref, - const char* rotate_query, const char* rotate_ref, - const bool run_Affine) -{ - - // Convert images to grayscale - Mat im1Gray, im2Gray; - cvtColor(im1, im1Gray, cv::COLOR_BGR2GRAY); - cvtColor(im2, im2Gray, cv::COLOR_BGR2GRAY); - - // Variables to store keypoints and descriptors - std::vector keypoints1, keypoints2; - Mat descriptors1, descriptors2; - - // Process images - Mat im1Proc, im2Proc, im1NormalProc; - im1Proc = preprocessImage(im1Gray, invert_query, flipflop_query, rotate_query); - im1NormalProc = preprocessImage(im1, FALSE, flipflop_query, rotate_query); - im2Proc = preprocessImage(im2Gray, invert_ref, flipflop_ref, rotate_ref); - - // Detect ORB features and compute descriptors. - Ptr orb = ORB::create(MAX_FEATURES); - orb->detectAndCompute(im1Proc, Mat(), keypoints1, descriptors1); - orb->detectAndCompute(im2Proc, Mat(), keypoints2, descriptors2); - Rcout << "DONE: orb based key-points detection and descriptors computation" << endl; - - // Match features. - std::vector matches; - Ptr matcher = DescriptorMatcher::create("BruteForce-Hamming"); - matcher->match(descriptors1, descriptors2, matches, Mat()); - Rcout << "DONE: BruteForce-Hamming - descriptor matching" << endl; - - // Sort matches by score - std::sort(matches.begin(), matches.end()); - - // Remove not so good matches - const int numGoodMatches = matches.size() * GOOD_MATCH_PERCENT; - matches.erase(matches.begin()+numGoodMatches, matches.end()); - Rcout << "DONE: get good matches by distance thresholding" << endl; - - // Extract location of good matches - std::vector points1, points2; - for( size_t i = 0; i < matches.size(); i++ ) - { - points1.push_back( keypoints1[ matches[i].queryIdx ].pt ); - points2.push_back( keypoints2[ matches[i].trainIdx ].pt ); - } - - // check variable - Rcout << "Calculating" << (run_Affine ? " (Affine) " : " (Homography) ") << "Transformation Matrix" << endl; - - // Find transformation matrix - cv::Mat mask; - if(run_Affine){ - std::vector match_mask; - h = estimateAffine2D(points1, - points2, - match_mask, - cv::RANSAC); - mask = IntVectorToMat(match_mask); - } else { - h = findHomography(points1, - points2, - cv::RANSAC, - 5, - mask); - } - - // Draw top matches and good ones only - std::vector top_matches; - std::vector keypoints1_best, keypoints2_best; - for(size_t i = 0; i < matches.size(); i++ ) - { - keypoints1_best.push_back(keypoints1[matches[i].queryIdx]); - keypoints2_best.push_back(keypoints2[matches[i].trainIdx]); - } - std::vector keypoints1_best2, keypoints2_best2; - int j=0; - for (int i = 0; i < mask.rows; i++) { - if (mask.at(i)) { - keypoints1_best2.push_back(keypoints1_best[i]); - keypoints2_best2.push_back(keypoints2_best[i]); - top_matches.push_back(cv::DMatch(static_cast(j), static_cast(j), 0)); - j++; - } - } - scaledDrawMatches(im1Proc, keypoints1_best2, im2Proc, keypoints2_best2, top_matches, imMatches); - - // Use homography to warp image - Mat im1Warp, im1NormalWarp; - if(h.rows == 2){ - warpAffine(im1Proc, im1Warp, h, im2Proc.size()); - warpAffine(im1NormalProc, im1NormalWarp, h, im2Proc.size()); - } else { - warpPerspective(im1Proc, im1Warp, h, im2Proc.size()); - warpPerspective(im1NormalProc, im1NormalWarp, h, im2Proc.size()); - } - - // Reverse process - im1Reg = reversepreprocessImage(im1NormalWarp, flipflop_ref, rotate_ref); - - // return as rgb - cvtColor(im2Proc, im2, cv::COLOR_GRAY2BGR); - - // resize image to visualize faster later in Shiny - im2 = resize_image(im2, 500); - im1Overlay = resize_image(im1Reg, 500); -} - -// align images with FLANN algorithm -void alignImagesFLANN(Mat &im1, Mat &im2, Mat &im1Reg, Mat &im1Overlay, - Mat &imMatches, Mat &h, Rcpp::List &keypoints, - const bool invert_query, const bool invert_ref, - const char* flipflop_query, const char* flipflop_ref, - const char* rotate_query, const char* rotate_ref, - const bool run_Affine, const bool run_TPS) -{ - - // parameters - cv::setRNGSeed(0); - SIFTParameters params; - - ////////////////////// - /// Process Images /// - ////////////////////// - - // Convert images to grayscale - Mat im1Gray, im2Gray; - cvtColor(im1, im1Gray, cv::COLOR_BGR2GRAY); - cvtColor(im2, im2Gray, cv::COLOR_BGR2GRAY); - - // Process images - Mat im1Proc, im2Proc, im1NormalProc; - im1Proc = preprocessImage(im1Gray, invert_query, flipflop_query, rotate_query); - im1NormalProc = preprocessImage(im1, FALSE, flipflop_query, rotate_query); - im2Proc = preprocessImage(im2Gray, invert_ref, flipflop_ref, rotate_ref); - - // //////////////////////////////////// - // /// Compute SIFT+FLANN+Homograpy /// - // //////////////////////////////////// - - // RUN SIFT+FLANN+Homography with retry - bool is_faulty = FALSE; - cv::Mat mask; - std::vector points1, points2; - // getSIFTTransformationMatrix(im1Proc, im2Proc, im1, im2, h, mask, imMatches, - // points1, points2, run_Affine, params, is_faulty); - - // check result - is_faulty = check_transformation_metrics(points1, points2, im2, h, mask); - Rcout << "MESSAGE: Registration is " << (is_faulty ? "degenerate!" : "not degenerate!") << endl; - - // Use homography to warp image - Mat im1Warp, im1NormalWarp; - if(h.rows == 2){ - warpAffine(im1Proc, im1Warp, h, im2Proc.size()); - warpAffine(im1NormalProc, im1NormalWarp, h, im2Proc.size()); - } else { - warpPerspective(im1Proc, im1Warp, h, im2Proc.size()); - warpPerspective(im1NormalProc, im1NormalWarp, h, im2Proc.size()); - } - - Rcout << "DONE: warped query image" << endl; - - /////////////////////// - /// Find Homography /// - /////////////////////// - - // continue with TPS or do FLANN only - Mat im1Reg_Warp_nonrigid; - Mat im1Reg_NormalWarp_nonrigid; - Mat im1Combine; - if(is_faulty || !run_TPS){ - - // change color map - cv::addWeighted(im2Proc, 0.7, im1Warp, 0.3, 0, im1Combine); - - // Reverse process - im1Reg = reversepreprocessImage(im1NormalWarp, flipflop_ref, rotate_ref); - - // return as rgb - cvtColor(im1Combine, im1Overlay, cv::COLOR_GRAY2BGR); - cvtColor(im2Proc, im2, cv::COLOR_GRAY2BGR); - - // TPS is requested (only if FLANN succeeded) - } else { - - Rcout << "MESSAGE: Running Thin-Plate-Spline Alignment" << endl; - - // Filtered points (inliers) based on the mask - std::vector filtered_points1; - std::vector filtered_points2; - for (int i = 0; i < mask.rows; i++) { - if (mask.at(i)) { - filtered_points1.push_back(points1[i]); - filtered_points2.push_back(points2[i]); - } - } - removeCloseMatches(filtered_points1, filtered_points2); - - // transform query - std::vector filtered_points1_reg; - if (h.rows == 2){ - cv::transform(filtered_points1, filtered_points1_reg, h); - } else { - cv::perspectiveTransform(filtered_points1, filtered_points1_reg, h); - } - - // get TPS matches - std::vector matches; - for (unsigned int i = 0; i < filtered_points2.size(); i++) - matches.push_back(cv::DMatch(i, i, 0)); - - // calculate TPS transformation - Ptr tps = cv::createThinPlateSplineShapeTransformer(0); - tps->estimateTransformation(filtered_points2, filtered_points1_reg, matches); - - // save keypoints - keypoints[0] = point2fToNumericMatrix(filtered_points2); - keypoints[1] = point2fToNumericMatrix(filtered_points1_reg); - - // determine extension limits for both images - int y_max = max(im1Warp.rows, im2.rows); - int x_max = max(im1Warp.cols, im2.cols); - - // extend images - cv::copyMakeBorder(im1Warp, im1Warp, 0.0, (int) (y_max - im1Warp.rows), 0.0, (x_max - im1Warp.cols), cv::BORDER_CONSTANT, Scalar(0, 0, 0)); - cv::copyMakeBorder(im1NormalWarp, im1NormalWarp, 0.0, (int) (y_max - im1NormalWarp.rows), 0.0, (x_max - im1NormalWarp.cols), cv::BORDER_CONSTANT, Scalar(0, 0, 0)); - - // transform image - Mat im1Reg_Warp_nonrigid; - Mat im1Reg_NormalWarp_nonrigid; - tps->warpImage(im1Warp, im1Reg_Warp_nonrigid); - tps->warpImage(im1NormalWarp, im1Reg_NormalWarp_nonrigid); - - // resize image - cv::Mat im1Reg_NormalWarp_nonrigid_cropped = im1Reg_NormalWarp_nonrigid(cv::Range(0,im2Proc.size().height), cv::Range(0,im2Proc.size().width)); - im1Reg_NormalWarp_nonrigid = im1Reg_NormalWarp_nonrigid_cropped.clone(); - - cv::Mat im1Reg_Warp_nonrigid_cropped = im1Reg_Warp_nonrigid(cv::Range(0,im2Proc.size().height), cv::Range(0,im2Proc.size().width)); - im1Reg_Warp_nonrigid = im1Reg_Warp_nonrigid_cropped.clone(); - - // change color map - cv::addWeighted(im2Proc, 0.7, im1Reg_Warp_nonrigid, 0.3, 0, im1Combine); - - // Reverse process - im1Reg = reversepreprocessImage(im1Reg_NormalWarp_nonrigid, flipflop_ref, rotate_ref); - - // return as rgb - cvtColor(im1Combine, im1Overlay, cv::COLOR_GRAY2BGR); - cvtColor(im2Proc, im2, cv::COLOR_GRAY2BGR); - } - - // resize image to visualize faster later in Shiny - im2 = resize_image(im2, 500); - im1Overlay = resize_image(im1Overlay, 500); -} - -// align images with FLANN algorithm -void alignImagesFLANN2(Mat &im1, Mat &im2, Mat &im1Reg, Mat &im1Overlay, - Mat &imMatches, Mat &h, Rcpp::List &keypoints, - const bool invert_query, const bool invert_ref, - const char* flipflop_query, const char* flipflop_ref, - const char* rotate_query, const char* rotate_ref, - const bool run_Affine, const bool run_TPS) -{ - - // parameters - cv::setRNGSeed(0); - SIFTParameters params; - - ////////////////////// - /// Process Images /// - ////////////////////// - - // Convert images to grayscale - Mat im1Gray, im2Gray; - cvtColor(im1, im1Gray, cv::COLOR_BGR2GRAY); - cvtColor(im2, im2Gray, cv::COLOR_BGR2GRAY); - - // Process images - Mat im1Proc, im2Proc, im1NormalProc; - im1Proc = preprocessImage(im1Gray, invert_query, flipflop_query, rotate_query); - im1NormalProc = preprocessImage(im1, FALSE, flipflop_query, rotate_query); - im2Proc = preprocessImage(im2Gray, invert_ref, flipflop_ref, rotate_ref); - - // //////////////////////////////////// - // /// Compute SIFT+FLANN+Homograpy /// - // //////////////////////////////////// - - // RUN SIFT+FLANN+Homography with retry - bool is_faulty = FALSE; - cv::Mat mask; - std::vector points1, points2; - - - // Variables to store keypoints and descriptors - std::vector keypoints1, keypoints2; - Mat descriptors1, descriptors2; - - // Detect SIFT features - // Ptr sift = cv::SIFT::create(params.sift_nfeatures); - Ptr sift = cv::SIFT::create(); - // computeSIFTTiles(im1Proc, keypoints1, descriptors1, sift, params); - // computeSIFTTiles(im2Proc, keypoints2, descriptors2, sift, params); - sift->detectAndCompute(im1Proc, Mat(), keypoints1, descriptors1); - sift->detectAndCompute(im2Proc, Mat(), keypoints2, descriptors2); - - Rcout << "MESSAGE: Generated " << keypoints1.size() << " and " << keypoints2.size() << " keypoints" << endl; - Rcout << "DONE: SIFT based key-points detection and descriptors computation" << endl; - - /////////////////////// - /// Compute FLANN ///// - /////////////////////// - - // Match features using FLANN matching - std::vector> matches; - cv::FlannBasedMatcher custom_matcher = cv::FlannBasedMatcher(cv::makePtr(5), cv::makePtr(50, 0, TRUE)); - cv::Ptr matcher = custom_matcher.create(); - matcher->knnMatch(descriptors1, descriptors2, matches, 2); - Rcout << "DONE: FLANN - Fast Library for Approximate Nearest Neighbors - descriptor matching" << endl; - - // Find good matches - // goodMatches = get_good_matches(matches) - std::vector good_matches; - getGoodMatches_temp(matches, good_matches); - Rcout << "DONE: get good matches by distance thresholding" << endl; - - /////////////////////// - /// Find Homography /// - /////////////////////// - - // Extract location of good matches - for( size_t i = 0; i < good_matches.size(); i++ ) - { - points1.push_back(keypoints1[good_matches[i].queryIdx].pt); - points2.push_back(keypoints2[good_matches[i].trainIdx].pt); - } - - // check variable - Rcout << "MESSAGE: Calculating" << (run_Affine ? " (Affine) " : " (Homography) ") << "Transformation Matrix" << endl; - - // Find transformation matrix - Rcout << "MESSAGE: Matching " << points1.size() << " keypoints" << endl; - if(run_Affine){ - std::vector match_mask; - h = estimateAffine2D(points1, - points2, - match_mask, - cv::RANSAC); - mask = IntVectorToMat(match_mask); - } else { - h = findHomography(points1, points2, RANSAC); - } - - // Draw top matches and good ones only - std::vector top_matches; - std::vector keypoints1_best, keypoints2_best; - for(size_t i = 0; i < good_matches.size(); i++ ) - { - keypoints1_best.push_back(keypoints1[good_matches[i].queryIdx]); - keypoints2_best.push_back(keypoints2[good_matches[i].trainIdx]); - } - std::vector keypoints1_best2, keypoints2_best2; - int j=0; - for (int i = 0; i < mask.rows; i++) { - if (mask.at(i)) { - keypoints1_best2.push_back(keypoints1_best[i]); - keypoints2_best2.push_back(keypoints2_best[i]); - top_matches.push_back(cv::DMatch(static_cast(j), static_cast(j), 0)); - j++; - } - } - // scaledDrawMatches(im1Proc, keypoints1_best2, im2Proc, keypoints2_best2, top_matches, imMatches); - drawMatches(im1Proc, keypoints1_best2, im2Proc, keypoints2_best2, top_matches, imMatches); - - // check result - is_faulty = check_transformation_metrics(points1, points2, im2, h, mask); - Rcout << "MESSAGE: Registration is " << (is_faulty ? "degenerate!" : "not degenerate!") << endl; - - // Use homography to warp image - Mat im1Warp, im1NormalWarp; - if(h.rows == 2){ - warpAffine(im1Proc, im1Warp, h, im2Proc.size()); - warpAffine(im1NormalProc, im1NormalWarp, h, im2Proc.size()); - } else { - warpPerspective(im1Proc, im1Warp, h, im2Proc.size()); - warpPerspective(im1NormalProc, im1NormalWarp, h, im2Proc.size()); - } - - Rcout << "DONE: warped query image" << endl; - - /////////////////////// - /// Find Homography /// - /////////////////////// - - // continue with TPS or do FLANN only - Mat im1Reg_Warp_nonrigid; - Mat im1Reg_NormalWarp_nonrigid; - Mat im1Combine; - if(is_faulty || !run_TPS){ - - // change color map - cv::addWeighted(im2Proc, 0.7, im1Warp, 0.3, 0, im1Combine); - - // Reverse process - im1Reg = reversepreprocessImage(im1NormalWarp, flipflop_ref, rotate_ref); - - // return as rgb - cvtColor(im1Combine, im1Overlay, cv::COLOR_GRAY2BGR); - cvtColor(im2Proc, im2, cv::COLOR_GRAY2BGR); - - // TPS is requested (only if FLANN succeeded) - } else { - - Rcout << "MESSAGE: Running Thin-Plate-Spline Alignment" << endl; - - // Filtered points (inliers) based on the mask - std::vector filtered_points1; - std::vector filtered_points2; - for (int i = 0; i < mask.rows; i++) { - if (mask.at(i)) { - filtered_points1.push_back(points1[i]); - filtered_points2.push_back(points2[i]); - } - } - removeCloseMatches(filtered_points1, filtered_points2); - - // transform query - std::vector filtered_points1_reg; - if (h.rows == 2){ - cv::transform(filtered_points1, filtered_points1_reg, h); - } else { - cv::perspectiveTransform(filtered_points1, filtered_points1_reg, h); - } - - // get TPS matches - std::vector matches; - for (unsigned int i = 0; i < filtered_points2.size(); i++) - matches.push_back(cv::DMatch(i, i, 0)); - - // calculate TPS transformation - Ptr tps = cv::createThinPlateSplineShapeTransformer(0); - tps->estimateTransformation(filtered_points2, filtered_points1_reg, matches); - - // save keypoints - keypoints[0] = point2fToNumericMatrix(filtered_points2); - keypoints[1] = point2fToNumericMatrix(filtered_points1_reg); - - // determine extension limits for both images - int y_max = max(im1Warp.rows, im2.rows); - int x_max = max(im1Warp.cols, im2.cols); - - // extend images - cv::copyMakeBorder(im1Warp, im1Warp, 0.0, (int) (y_max - im1Warp.rows), 0.0, (x_max - im1Warp.cols), cv::BORDER_CONSTANT, Scalar(0, 0, 0)); - cv::copyMakeBorder(im1NormalWarp, im1NormalWarp, 0.0, (int) (y_max - im1NormalWarp.rows), 0.0, (x_max - im1NormalWarp.cols), cv::BORDER_CONSTANT, Scalar(0, 0, 0)); - - // transform image - Mat im1Reg_Warp_nonrigid; - Mat im1Reg_NormalWarp_nonrigid; - tps->warpImage(im1Warp, im1Reg_Warp_nonrigid); - tps->warpImage(im1NormalWarp, im1Reg_NormalWarp_nonrigid); - - // resize image - cv::Mat im1Reg_NormalWarp_nonrigid_cropped = im1Reg_NormalWarp_nonrigid(cv::Range(0,im2Proc.size().height), cv::Range(0,im2Proc.size().width)); - im1Reg_NormalWarp_nonrigid = im1Reg_NormalWarp_nonrigid_cropped.clone(); - - cv::Mat im1Reg_Warp_nonrigid_cropped = im1Reg_Warp_nonrigid(cv::Range(0,im2Proc.size().height), cv::Range(0,im2Proc.size().width)); - im1Reg_Warp_nonrigid = im1Reg_Warp_nonrigid_cropped.clone(); - - // change color map - cv::addWeighted(im2Proc, 0.7, im1Reg_Warp_nonrigid, 0.3, 0, im1Combine); - - // Reverse process - im1Reg = reversepreprocessImage(im1Reg_NormalWarp_nonrigid, flipflop_ref, rotate_ref); - - // return as rgb - cvtColor(im1Combine, im1Overlay, cv::COLOR_GRAY2BGR); - cvtColor(im2Proc, im2, cv::COLOR_GRAY2BGR); - } - - // resize image to visualize faster later in Shiny - im2 = resize_image(im2, 500); - im1Overlay = resize_image(im1Overlay, 500); } \ No newline at end of file diff --git a/src/auxiliary.cpp b/src/auxiliary.cpp index ba4a455f..4a128a35 100644 --- a/src/auxiliary.cpp +++ b/src/auxiliary.cpp @@ -35,81 +35,12 @@ Rcpp::NumericMatrix replaceNaMatrix(Rcpp::NumericMatrix mat, int replace) { return mat; } -//// -// memory -//// - -// // memory check -// void log_mem_usage(const std::string& label = "") { -// struct rusage usage; -// getrusage(RUSAGE_SELF, &usage); -// long rss_b = usage.ru_maxrss; -// -// double rss_kb = rss_b / 1024.0; -// double rss_mb = rss_kb / 1024.0; -// double rss_gb = rss_mb / 1024.0; -// -// Rcpp::Rcout << "Used Memory [" << label << "]: " << rss_gb << " GB" << std::endl; -// } -// -// void log_mem_macos(const std::string& label = "") { -// mach_task_basic_info info; -// mach_msg_type_number_t size = MACH_TASK_BASIC_INFO_COUNT; -// kern_return_t kr = task_info(mach_task_self(), MACH_TASK_BASIC_INFO, -// (task_info_t)&info, &size); -// -// if (kr != KERN_SUCCESS) { -// Rcpp::Rcerr << "[MEM " << label << "] Failed to get memory info.\n"; -// return; -// } -// -// double rss_gb = static_cast(info.resident_size) / (1024.0 * 1024.0 * 1024.0); -// double virt_gb = static_cast(info.virtual_size) / (1024.0 * 1024.0 * 1024.0); -// -// Rcpp::Rcout << "[MEM " << label << "] Resident (RSS): " -// << rss_gb << " GB, Virtual: " << virt_gb << " GB\n"; -// } -// -// double object_size_long(long bsize) { -// -// double rss_kb = bsize / 1024.0; -// double rss_mb = rss_kb / 1024.0; -// double rss_gb = rss_mb / 1024.0; -// -// return rss_gb; -// } -// -// double object_size_double(double bsize) { -// -// double rss_kb = bsize / 1024; -// double rss_mb = rss_kb / 1024; -// double rss_gb = rss_mb / 1024; -// -// return rss_gb; -// } -// -// double get_resident_bytes() { -// mach_task_basic_info info; -// mach_msg_type_number_t size = MACH_TASK_BASIC_INFO_COUNT; -// if (task_info(mach_task_self(), MACH_TASK_BASIC_INFO, -// (task_info_t)&info, &size) != KERN_SUCCESS) { -// return 0; -// } -// return static_cast(info.resident_size); -// } -// -// double bytes_to_gb(double bytes) { -// return bytes / (1024.0 * 1024.0 * 1024.0); -// } - //// // Conversion //// // Function to convert a cv::Mat object to a RawVector for magick images Rcpp::RawVector matToImage(const cv::Mat &mat) { - // profiler - // MemProfiler mp("Mat -> Image"); // Create RawVector object Rcpp::RawVector rawvec(mat.total() * mat.elemSize()); @@ -123,8 +54,6 @@ Rcpp::RawVector matToImage(const cv::Mat &mat) { // Function to convert a RawVector for magick images to a cv::Mat object cv::Mat imageToMat(Rcpp::RawVector &image_data, int width, int height) { - // profiler - // MemProfiler mp("Image -> Mat"); // Create cv::Mat object cv::Mat mat(height, width, CV_8UC3, image_data.begin()); @@ -135,6 +64,22 @@ cv::Mat imageToMat(Rcpp::RawVector &image_data, int width, int height) { return mat; } +// Function to convert a cv::Mat object to a RawVector for magick images +Rcpp::IntegerVector matToMask(const cv::Mat &mat) { + + cv::Mat intMat; + mat.convertTo(intMat, CV_32S); + Rcpp::IntegerVector intvec(intMat.total()); + std::memcpy( + intvec.begin(), + intMat.data, + static_cast(intvec.size()) * sizeof(int) + ); + intvec.attr("dim") = Rcpp::Dimension(intMat.rows, intMat.cols); + + return intvec; +} + // Function to convert a NumericMatrix object to a cv::Mat cv::Mat numericMatrixToMat(Rcpp::NumericMatrix nm) { cv::Mat m(nm.rows(), nm.cols(), CV_64F); @@ -197,6 +142,17 @@ std::vector Point2fToDoubleVector(std::vector &points) { return vec; } +// Function to convert a cv::Keypoint object to a std::vector +std::vector KeyPointToPoint2f(std::vector &keypoints) { + int n = keypoints.size(); + std::vector points; + + for (int i = 0; i < n; i++) { + points.push_back(keypoints[i].pt); + } + return points; +} + // Function to convert a cv::Point2f object to a cv::Mat std::vector matToPoint2f(cv::Mat &mat) { std::vector points; @@ -289,3 +245,61 @@ double cppSD(std::vector &points) std::vector().swap(inVec); return std::sqrt( sd / (n-1) ); } + +double meanDistances(std::vector& pts1, + std::vector& pts2) +{ + if (pts1.size() != pts2.size() || pts1.empty()) + return 0.0; + + double sumDist = 0.0; + for (size_t i = 0; i < pts1.size(); ++i) + { + const double dx = pts1[i].x - pts2[i].x; + const double dy = pts1[i].y - pts2[i].y; + sumDist += std::sqrt(dx * dx + dy * dy); + } + + return sumDist / pts1.size(); +} + +double medianDistances(std::vector& pts1, + std::vector& pts2) +{ + if (pts1.size() != pts2.size() || pts1.empty()) + return 0.0; + + std::vector distances; + distances.reserve(pts1.size()); + + for (size_t i = 0; i < pts1.size(); ++i) + { + const double dx = pts1[i].x - pts2[i].x; + const double dy = pts1[i].y - pts2[i].y; + distances.push_back(std::sqrt(dx * dx + dy * dy)); + } + + const size_t n = distances.size(); + const size_t mid = n / 2; + + std::nth_element(distances.begin(), + distances.begin() + mid, + distances.end()); + + if (n % 2 == 1) + { + return distances[mid]; + } + else + { + double upper = distances[mid]; + + std::nth_element(distances.begin(), + distances.begin() + mid - 1, + distances.end()); + + double lower = distances[mid - 1]; + + return (lower + upper) / 2.0; + } +} \ No newline at end of file diff --git a/src/auxiliary.h b/src/auxiliary.h index 0985e819..9d9d5271 100644 --- a/src/auxiliary.h +++ b/src/auxiliary.h @@ -19,6 +19,10 @@ Rcpp::NumericMatrix replaceNaMatrix(Rcpp::NumericMatrix mat, int replace); Rcpp::RawVector matToImage(const cv::Mat &mat); cv::Mat imageToMat(Rcpp::RawVector &image_data, int width, int height); +// cv::Mat vs Rcpp::RawVector(Image) with 2 dim (mostly for masks) +Rcpp::IntegerVector matToMask(const cv::Mat &mat); +// Rcpp::RawVector matToMask(const cv::Mat &mat); + // cv::Mat vs Rcpp::NumericMatrix cv::Mat numericMatrixToMat(Rcpp::NumericMatrix nm); Rcpp::NumericMatrix matToNumericMatrix(cv::Mat m); @@ -37,6 +41,9 @@ cv::Mat IntVectorToMat(std::vector &points); // std::vector vs std::vector std::vector KeyPointToDoubleVector(std::vector &points); std::vector Point2fToDoubleVector(std::vector &points); + +// std::vector vs std::vector +std::vector KeyPointToPoint2f(std::vector &keypoints); //// // stats @@ -46,32 +53,8 @@ std::vector Point2fToDoubleVector(std::vector &points); double cppSD(std::vector &points); double cppSD(std::vector &points); -//// -// memory -//// - -// void log_mem_usage(const std::string& label); -// void log_mem_macos(const std::string& label); -// double object_size_long(long bsize); -// double object_size_double(double bsize); -// double get_resident_bytes(); -// double bytes_to_gb(double bytes); -// -// struct MemProfiler { -// size_t start; -// std::string label; -// -// MemProfiler(const std::string& lbl) : label(lbl) { -// start = get_resident_bytes(); -// } -// -// ~MemProfiler() { -// size_t end = get_resident_bytes(); -// double diff = (double) end - (double) start; -// if(diff < 0.0) diff = 0.0; -// double diff_gb = bytes_to_gb(diff); -// Rcpp::Rcout << "[MEM] " << label << ": +" << diff_gb << " GB" << std::endl; -// } -// }; +// mean distance between points +double meanDistances(std::vector &pts1, std::vector &pts2); +double medianDistances(std::vector &pts1, std::vector &pts2); #endif \ No newline at end of file diff --git a/src/image.cpp b/src/image.cpp index 5152d859..030b23a8 100644 --- a/src/image.cpp +++ b/src/image.cpp @@ -126,18 +126,41 @@ void scaledDrawMatches(cv::Mat im1, std::vector &keypoints1, // draw matches drawMatches(im1, keypoints1, im2, keypoints2, top_matches, imMatches); } + +cv::Mat warpTPSImage(cv::Mat& ref_image, + cv::Mat& query_image, + Ptr& tps, + const int border_y, + const int border_x, + const int interpolation){ -// [[Rcpp::export]] -Rcpp::RawVector warpRcppImage(Rcpp::RawVector ref_image, Rcpp::RawVector query_image, - Rcpp::List mapping, - const int width1, const int height1, - const int width2, const int height2) -{ - // Read reference image - cv::Mat imReference = imageToMat(ref_image, width1, height1); + // determine extension limits for both images + int y_max = max(query_image.rows, border_y); + int x_max = max(query_image.cols, border_x); + + // extend images + cv::copyMakeBorder(query_image, query_image, + 0.0, (int) (y_max - query_image.rows), + 0.0, (x_max - query_image.cols), + cv::BORDER_CONSTANT, Scalar(0, 0, 0)); + + // transform image + cv::Mat query_image_reg; + tps->warpImage(query_image, query_image_reg, interpolation); + + // resize image + query_image_reg = query_image_reg( + cv::Range(0,ref_image.size().height), + cv::Range(0,ref_image.size().width)); + + // return + return query_image_reg.clone(); +} + +void warpImage(cv::Mat& ref_image, + cv::Mat& query_image, + Rcpp::List mapping){ - // Read image to be aligned - cv::Mat im = imageToMat(query_image, width2, height2); cv::Mat im_temp; // list @@ -155,11 +178,11 @@ Rcpp::RawVector warpRcppImage(Rcpp::RawVector ref_image, Rcpp::RawVector query_i // transform coordinates if(h.rows == 2){ - cv::warpAffine(im, im_temp, h, imReference.size()); + cv::warpAffine(query_image, im_temp, h, ref_image.size()); } else { - cv::warpPerspective(im, im_temp, h, imReference.size()); + cv::warpPerspective(query_image, im_temp, h, ref_image.size()); } - im = im_temp; + query_image = im_temp; } // non-rigid warping @@ -177,35 +200,64 @@ Rcpp::RawVector warpRcppImage(Rcpp::RawVector ref_image, Rcpp::RawVector query_i // calculate transformation Ptr tps = cv::createThinPlateSplineShapeTransformer(0); - tps->estimateTransformation(ref_mat, query_mat, matches); - // determine extension limits for both images - int y_max = max(im.rows, imReference.rows); - int x_max = max(im.cols, imReference.cols); - - // extend images - cv::copyMakeBorder(im, im, 0.0, (int) (y_max - im.rows), 0.0, (x_max - im.cols), cv::BORDER_CONSTANT, Scalar(0, 0, 0)); - - // transform image - tps->warpImage(im, im_temp); + // estimate transformation + tps->estimateTransformation(ref_mat, query_mat, matches); - // resize image - cv::Mat im_temp_cropped = im_temp(cv::Range(0,imReference.size().height), cv::Range(0,imReference.size().width)); - im_temp = im_temp_cropped.clone(); + // transform image using trained tps + im_temp = warpTPSImage(ref_image, query_image, tps, + ref_image.rows, ref_image.cols, cv::INTER_LINEAR); + + // // determine extension limits for both images + // int y_max = max(query_image.rows, ref_image.rows); + // int x_max = max(query_image.cols, ref_image.cols); + // + // // extend images + // cv::copyMakeBorder(query_image, query_image, + // 0.0, (int) (y_max - query_image.rows), + // 0.0, (x_max - query_image.cols), + // cv::BORDER_CONSTANT, Scalar(0, 0, 0)); + // + // // transform image + // tps->warpImage(query_image, im_temp); + // + // // resize image + // cv::Mat im_temp_cropped = im_temp(cv::Range(0,ref_image.size().height), cv::Range(0,ref_image.size().width)); + // im_temp = im_temp_cropped.clone(); } else { // pass registered object - im_temp = im; + im_temp = query_image; } - im = im_temp; + query_image = im_temp; } +} + +// [[Rcpp::export]] +Rcpp::RawVector warpRcppImage(Rcpp::RawVector ref_image, + Rcpp::RawVector query_image, + Rcpp::List mapping, + const int width1, const int height1, + const int width2, const int height2) +{ + // Read reference image + cv::Mat imReference = imageToMat(ref_image, width1, height1); + + // Read image to be aligned + cv::Mat im = imageToMat(query_image, width2, height2); + + // warp image + warpImage(imReference, im, mapping); // return return matToImage(im); } +///// +// Legacy //// +///// // [[Rcpp::export]] Rcpp::RawVector warpImageAuto(Rcpp::RawVector ref_image, Rcpp::RawVector query_image, diff --git a/src/image.h b/src/image.h index 6a9a7b26..34fe97ef 100644 --- a/src/image.h +++ b/src/image.h @@ -2,6 +2,10 @@ #include #include "opencv2/shape/shape_transformer.hpp" +using namespace Rcpp; +using namespace std; +using namespace cv; + #ifndef IMAGE_H #define IMAGE_H @@ -9,8 +13,10 @@ // Processing //// -cv::Mat preprocessImage(cv::Mat &im, const bool invert, const char* flipflop, const char* rotate); -cv::Mat reversepreprocessImage(cv::Mat &im, const char* flipflop, const char* rotate); +cv::Mat preprocessImage(cv::Mat &im, const bool invert, + const char* flipflop, const char* rotate); +cv::Mat reversepreprocessImage(cv::Mat &im, + const char* flipflop, const char* rotate); cv::Mat resize_image(cv::Mat &im, int width); std::vector resize_keypoints(std::vector &keypoints, cv::Mat &im, @@ -24,17 +30,27 @@ void scaledDrawMatches(cv::Mat im1, std::vector &keypoints1, // Warping //// -Rcpp::RawVector warpImage(Rcpp::RawVector ref_image, Rcpp::RawVector query_image, +cv::Mat warpTPSImage(cv::Mat& ref_image, + cv::Mat& query_image, + Ptr& tps, + const int border_x, + const int border_y, + const int interpolation); + +Rcpp::RawVector warpImage(Rcpp::RawVector ref_image, + Rcpp::RawVector query_image, Rcpp::List mapping, const int width1, const int height1, const int width2, const int height2); -Rcpp::RawVector warpImageAuto(Rcpp::RawVector ref_image, Rcpp::RawVector query_image, +Rcpp::RawVector warpImageAuto(Rcpp::RawVector ref_image, + Rcpp::RawVector query_image, Rcpp::List mapping, const int width1, const int height1, const int width2, const int height2); -Rcpp::RawVector warpImageManual(Rcpp::RawVector ref_image, Rcpp::RawVector query_image, +Rcpp::RawVector warpImageManual(Rcpp::RawVector ref_image, + Rcpp::RawVector query_image, Rcpp::List mapping, const int width1, const int height1, const int width2, const int height2); diff --git a/src/manual_registration.cpp b/src/manual_registration.cpp index 7cd978bb..1474ee82 100644 --- a/src/manual_registration.cpp +++ b/src/manual_registration.cpp @@ -4,8 +4,11 @@ #include #include "opencv2/shape/shape_transformer.hpp" -// Auxiliary +// Library #include "auxiliary.h" +#include "image.h" +#include "metrics.h" +#include "matte_mi.h" // Namespaces using namespace Rcpp; @@ -14,7 +17,10 @@ using namespace cv; // align images with TPS algorithm void alignImagesTPS(Mat &im1, Mat &im2, Mat &im1Reg, Rcpp::List &keypoints, - Rcpp::NumericMatrix query_landmark, Rcpp::NumericMatrix reference_landmark) + Rcpp::NumericMatrix query_landmark, Rcpp::NumericMatrix reference_landmark, + const bool invert_query, const bool invert_ref, + Mat1d &accuracyMatte, + std::map &accuracy) { // seed @@ -31,6 +37,9 @@ void alignImagesTPS(Mat &im1, Mat &im2, Mat &im1Reg, Rcpp::List &keypoints, for (unsigned int i = 0; i < ref_mat.size(); i++) matches.push_back(cv::DMatch(i, i, 0)); + // message + Rcout << "Running Coarse Alignment (Thin-Plate-Spline)" << endl; + // calculate transformation Ptr tps = cv::createThinPlateSplineShapeTransformer(0); tps->estimateTransformation(ref_mat, query_mat, matches); @@ -39,19 +48,25 @@ void alignImagesTPS(Mat &im1, Mat &im2, Mat &im1Reg, Rcpp::List &keypoints, keypoints[0] = point2fToNumericMatrix(ref_mat); keypoints[1] = point2fToNumericMatrix(query_mat); - // determine extension limits for both images - int y_max = max(im1.rows, im2.rows); - int x_max = max(im1.cols, im2.cols); - - // extend images - cv::copyMakeBorder(im1, im1, 0.0, (int) (y_max - im1.rows), 0.0, (x_max - im1.cols), cv::BORDER_CONSTANT, Scalar(0, 0, 0)); - - // transform image - tps->warpImage(im1, im1Reg); - - // resize image - cv::Mat im1Reg_cropped = im1Reg(cv::Range(0,im2.size().height), cv::Range(0,im2.size().width)); - im1Reg = im1Reg_cropped.clone(); + // transform image using trained tps + im1Reg = warpTPSImage(im2, im1, tps, + im2.rows, im2.cols, cv::INTER_LINEAR); + + // process + Mat im1Proc, im2Proc; + cvtColor(im1Reg, im1Proc, cv::COLOR_BGR2GRAY); + cvtColor(im2, im2Proc, cv::COLOR_BGR2GRAY); + im1Proc = preprocessImage(im1Proc, invert_query, "None", "0"); + im2Proc = preprocessImage(im2Proc, invert_ref, "None", "0"); + + // get alignment mask + cv::Mat alignmentMask = generateOverlapMask(im2Proc, + tps, + im1.size()); + + // get alignment metrics + accuracy = getAlignmentMetrics(im1Proc, im2Proc, alignmentMask, "Coarse"); + accuracyMatte = MatteMIMap(im2Proc, im1Proc, alignmentMask, 50); } // align images with TPS algorithm @@ -95,9 +110,21 @@ void alignImagesTPS_points(Rcpp::NumericMatrix &query_data, } // align images with FLANN algorithm -void alignImagesAffineTPS(Mat &im1, Mat &im2, Mat &im1Reg, Mat &h, Rcpp::List &keypoints, - Rcpp::NumericMatrix query_landmark, Rcpp::NumericMatrix reference_landmark, - const bool run_Affine, const bool run_TPS) +void alignImagesAffineTPS(Mat &im1, + Mat &im2, + Mat &im1Reg, + Mat &h, + Rcpp::List &keypoints, + Rcpp::NumericMatrix query_landmark, + Rcpp::NumericMatrix reference_landmark, + const bool invert_query, + const bool invert_ref, + const bool run_Affine, + const bool run_TPS, + const bool compute_matte_map, + Mat1d &accuracyMatte, + std::map &accuracy_coarse, + std::map &accuracy_fine) { // seed cv::setRNGSeed(0); @@ -115,6 +142,8 @@ void alignImagesAffineTPS(Mat &im1, Mat &im2, Mat &im1Reg, Mat &h, Rcpp::List &k // calculate homography transformation Rcout << "Calculating" << (run_Affine ? " (Affine) " : " (Homography) ") << "Transformation Matrix" << endl; + + // warp image Mat im1Affine; std::vector query_reg; if(run_Affine){ @@ -127,6 +156,20 @@ void alignImagesAffineTPS(Mat &im1, Mat &im2, Mat &im1Reg, Mat &h, Rcpp::List &k cv::perspectiveTransform(query_mat, query_reg, h); } + // get alignment metrics for Coarse registration + cv::Mat alignmentMask = generateOverlapMask(im2.size(), + h, + im1.size()); + + // get matte metric, process image before + Mat im1Proc, im2Proc; + cvtColor(im1Affine, im1Proc, cv::COLOR_BGR2GRAY); + cvtColor(im2, im2Proc, cv::COLOR_BGR2GRAY); + im1Proc = preprocessImage(im1Proc, invert_query, "None", "0"); + im2Proc = preprocessImage(im2Proc, invert_ref, "None", "0"); + accuracy_coarse = getAlignmentMetrics(im1Proc, im2Proc, alignmentMask, "Coarse"); + accuracyMatte = MatteMIMap(im2Proc, im1Proc, alignmentMask, 50); + if(!run_TPS){ // clone and exit @@ -135,8 +178,8 @@ void alignImagesAffineTPS(Mat &im1, Mat &im2, Mat &im1Reg, Mat &h, Rcpp::List &k } else { // message - Rcout << "Running Thin-Plate-Spline Alignment" << endl; - + Rcout << "Running Fine Alignment (Thin-Plate-Spline)" << endl; + // calculate TPS transformation Ptr tps = cv::createThinPlateSplineShapeTransformer(0); tps->estimateTransformation(ref_mat, query_reg, matches); @@ -145,19 +188,23 @@ void alignImagesAffineTPS(Mat &im1, Mat &im2, Mat &im1Reg, Mat &h, Rcpp::List &k keypoints[0] = point2fToNumericMatrix(ref_mat); keypoints[1] = point2fToNumericMatrix(query_reg); - // determine extension limits for both images - int y_max = max(im1Affine.rows, im2.rows); - int x_max = max(im1Affine.cols, im2.cols); + // warp overlap mask + alignmentMask = warpTPSImage(im2, alignmentMask, tps, + im2.rows, im2.cols, + cv::INTER_NEAREST); - // extend images - cv::copyMakeBorder(im1Affine, im1Affine, 0.0, (int) (y_max - im1Affine.rows), 0.0, (x_max - im1Affine.cols), cv::BORDER_CONSTANT, Scalar(0, 0, 0)); + // transform image using trained tps + im1Reg = warpTPSImage(im2, im1Affine, tps, + im2.rows, im2.cols, cv::INTER_LINEAR); - // transform image - tps->warpImage(im1Affine, im1Reg); - - // resize image - cv::Mat im1Reg_cropped = im1Reg(cv::Range(0,im2.size().height), cv::Range(0,im2.size().width)); - im1Reg = im1Reg_cropped.clone(); + // get matte metric, process + // im2 is already processed + Mat im1Proc; + cvtColor(im1Reg, im1Proc, cv::COLOR_BGR2GRAY); + im1Proc = preprocessImage(im1Proc, invert_query, "None", "0"); + accuracy_fine = getAlignmentMetrics(im1Proc, im2Proc, alignmentMask, "Fine"); + if(compute_matte_map) + accuracyMatte = MatteMIMap(im2Proc, im1Proc, alignmentMask, 50); } } @@ -175,6 +222,9 @@ void alignImagesAffineTPS_points(Rcpp::NumericMatrix &query_data, RNG rng(12345); Scalar value; + // message + Rcout << "Running Coarse Alignment (Manual)" << endl; + // Get landmarks as Point2f std::vector query_mat = numericMatrixToPoint2f(query_landmark); std::vector ref_mat = numericMatrixToPoint2f(reference_landmark); @@ -189,6 +239,7 @@ void alignImagesAffineTPS_points(Rcpp::NumericMatrix &query_data, // calculate homography transformation Rcout << "Calculating" << (run_Affine ? " (Affine) " : " (Homography) ") << "Transformation Matrix" << endl; + std::vector query_reg; std::vector query_data_reg; if(run_Affine){ @@ -204,8 +255,8 @@ void alignImagesAffineTPS_points(Rcpp::NumericMatrix &query_data, if(run_TPS){ // message - Rcout << "Running Thin-Plate-Spline Alignment" << endl; - + Rcout << "Running Fine Alignment (Thin-Plate-Spline)" << endl; + // calculate TPS transformation Ptr tps = cv::createThinPlateSplineShapeTransformer(0); tps->estimateTransformation(ref_mat, query_reg, matches); @@ -232,14 +283,19 @@ Rcpp::List manual_registeration_rawvector(Rcpp::RawVector ref_image, const int height1, const int width2, const int height2, + const bool invert_query, + const bool invert_ref, Rcpp::String method, - Rcpp::String nonrigid) + Rcpp::String nonrigid, + const bool compute_matte_map = true) { // Return data - Rcpp::List out(2); + Rcpp::List out(5); Rcpp::List out_trans(2); Rcpp::List keypoints(2); Mat imReg, h; + Mat1d accuracyMatte; + std::map accuracy_coarse, accuracy_fine; // get params const bool run_TPS = (strcmp(method.get_cstring(), "Homography + Non-Rigid") == 0 || @@ -257,14 +313,26 @@ Rcpp::List manual_registeration_rawvector(Rcpp::RawVector ref_image, alignImagesAffineTPS(im, imReference, imReg, h, keypoints, query_landmark, reference_landmark, - run_Affine, run_TPS); + invert_query, + invert_ref, + run_Affine, + run_TPS, + compute_matte_map, + accuracyMatte, + accuracy_coarse, + accuracy_fine); } // Non-rigid (TPS) only if(strcmp(method.get_cstring(), "Non-Rigid") == 0){ alignImagesTPS(im, imReference, imReg, keypoints, - query_landmark, reference_landmark); + query_landmark, + reference_landmark, + invert_query, + invert_ref, + accuracyMatte, + accuracy_coarse); } // transformation matrix, can be either a matrix, set of keypoints or both @@ -272,8 +340,11 @@ Rcpp::List manual_registeration_rawvector(Rcpp::RawVector ref_image, out_trans[1] = keypoints; out[0] = out_trans; - // registered image if exists + // registered image and accuracy if exists out[1] = matToImage(imReg.clone()); + out[2] = matToNumericMatrix(accuracyMatte); // Matte MI metric + out[3] = accuracy_coarse; + out[4] = accuracy_fine; return out; } @@ -313,8 +384,6 @@ Rcpp::List manual_registeration_matrix(Rcpp::NumericMatrix query_data, keypoints, query_landmark, reference_landmark); - keypoints[0] = keypoints[0]; - keypoints[1] = keypoints[1]; } // transformation matrix, can be either a matrix, set of keypoints or both diff --git a/src/mapping.cpp b/src/mapping.cpp index 2a10b90a..c3bbca0b 100644 --- a/src/mapping.cpp +++ b/src/mapping.cpp @@ -66,13 +66,6 @@ Rcpp::NumericMatrix applyRcppMapping(Rcpp::NumericMatrix coords, Rcpp::List mapp // apply transformation to coordinates tps->applyTransformation(coords_mat, coords_temp); - // // temp position - // std::vector query_mat_tmp; - // tps->applyTransformation(query_mat, query_mat_tmp); - // Rcout << query_mat << endl; - // Rcout << query_mat_tmp << endl; - // Rcout << ref_mat << endl; - } else { coords_temp = coords_mat; } diff --git a/src/matte_mi.cpp b/src/matte_mi.cpp new file mode 100644 index 00000000..aa55ba23 --- /dev/null +++ b/src/matte_mi.cpp @@ -0,0 +1,443 @@ +#include "matte_mi.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace { + +using Pixel = unsigned char; + +struct IntensityRange { + double min; + double max; +}; + +struct ChunkSize { + int height = 50; + int width = 50; +}; + +struct GlobalCounts { + std::array fixed{}; + std::array moving{}; + std::size_t validPairs = 0U; +}; + +// Source pixels remain CV_8U. All Mattes arithmetic remains double. +double cubicBSpline(double u) noexcept { + u = std::abs(u); + if (u < 1.0) { + const double u2 = u * u; + return (4.0 - 6.0 * u2 + 3.0 * u2 * u) / 6.0; + } + if (u < 2.0) { + const double t = 2.0 - u; + return t * t * t / 6.0; + } + return 0.0; +} + +double scaleToBinPosition( + double value, + IntensityRange range, + double low, + double high) { + if (!std::isfinite(range.min) || + !std::isfinite(range.max) || + !(range.max > range.min)) { + throw std::invalid_argument("Invalid intensity range."); + } + value = std::clamp(value, range.min, range.max); + return low + (value - range.min) * (high - low) / + (range.max - range.min); +} + +std::size_t roundToNearestEvenNonnegative(double x) noexcept { + const double lowerDouble = std::floor(x); + const double fraction = x - lowerDouble; + const auto lower = static_cast(lowerDouble); + if (fraction < 0.5) return lower; + if (fraction > 0.5) return lower + 1U; + return (lower % 2U == 0U) ? lower : lower + 1U; +} + +bool isValidRange(IntensityRange range) noexcept { + return std::isfinite(range.min) && + std::isfinite(range.max) && + range.max > range.min; +} + +int ceilDividePositive(int value, int divisor) noexcept { + return value / divisor + ((value % divisor) != 0 ? 1 : 0); +} + +void validateInputs( + const cv::Mat& fixed, + const cv::Mat& moving, + const cv::Mat& mask, + int bins) { + if (fixed.empty() || moving.empty()) { + throw std::invalid_argument( + "fixed and moving images must not be empty."); + } + if (fixed.type() != CV_8UC1 || moving.type() != CV_8UC1) { + throw std::invalid_argument( + "The zero-conversion implementation expects fixed and moving " + "to be CV_8UC1."); + } + if (fixed.size() != moving.size()) { + throw std::invalid_argument( + "fixed and moving must have the same dimensions."); + } + if (!mask.empty() && + (mask.type() != CV_8UC1 || mask.size() != fixed.size())) { + throw std::invalid_argument( + "mask must be empty or CV_8UC1 with the same dimensions."); + } + if (bins < 4) { + throw std::invalid_argument( + "bins must be >= 4 for cubic B-spline smoothing."); + } +} + +GlobalCounts collectGlobalCounts( + const cv::Mat& fixed, + const cv::Mat& moving, + const cv::Mat& mask) { + GlobalCounts out; + for (int y = 0; y < fixed.rows; ++y) { + const Pixel* fixedRow = fixed.ptr(y); + const Pixel* movingRow = moving.ptr(y); + const Pixel* maskRow = mask.empty() ? nullptr : mask.ptr(y); + for (int x = 0; x < fixed.cols; ++x) { + if (maskRow != nullptr && maskRow[x] == 0U) continue; + ++out.fixed[fixedRow[x]]; + ++out.moving[movingRow[x]]; + ++out.validPairs; + } + } + return out; +} + +Pixel valueAtRank( + const std::array& counts, + std::size_t rank) { + std::size_t cumulative = 0U; + for (std::size_t value = 0U; value < counts.size(); ++value) { + cumulative += counts[value]; + if (rank < cumulative) return static_cast(value); + } + throw std::out_of_range("Percentile rank is out of range."); +} + +// Exact NumPy-style linear percentile for CV_8U values, without sorting pixels. +double percentileFromCounts( + const std::array& counts, + std::size_t count, + double percentile) { + if (count == 0U) { + throw std::invalid_argument("Cannot calculate an empty percentile."); + } + if (!std::isfinite(percentile) || percentile < 0.0 || percentile > 100.0) { + throw std::invalid_argument("Percentile must be in [0, 100]."); + } + if (count == 1U) return static_cast(valueAtRank(counts, 0U)); + + const double index = static_cast(count - 1U) * percentile / 100.0; + const auto lowerIndex = static_cast(std::floor(index)); + const auto upperIndex = static_cast(std::ceil(index)); + const double fraction = index - static_cast(lowerIndex); + const double lower = static_cast(valueAtRank(counts, lowerIndex)); + const double upper = static_cast(valueAtRank(counts, upperIndex)); + return lower + fraction * (upper - lower); +} + +IntensityRange minMaxRangeFromCounts( + const std::array& counts, + std::size_t count) { + if (count == 0U) { + const double nan = std::numeric_limits::quiet_NaN(); + return {nan, nan}; + } + std::size_t minimum = 0U; + while (minimum < counts.size() && counts[minimum] == 0U) ++minimum; + std::size_t maximum = counts.size() - 1U; + while (maximum > 0U && counts[maximum] == 0U) --maximum; + return { + static_cast(minimum), + static_cast(maximum) + }; +} + +void addMattesPair( + double fixedValue, + double movingValue, + IntensityRange fixedRange, + IntensityRange movingRange, + std::size_t bins, + std::vector& jointHistogram) { + const double fixedPosition = scaleToBinPosition( + fixedValue, fixedRange, 0.0, static_cast(bins - 1U)); + std::size_t fixedBin = roundToNearestEvenNonnegative(fixedPosition); + fixedBin = std::min(fixedBin, bins - 1U); + + const double movingPosition = scaleToBinPosition( + movingValue, movingRange, 1.0, static_cast(bins - 2U)); + const auto baseBin = static_cast(std::floor(movingPosition)); + + for (int offset = -1; offset <= 2; ++offset) { + const std::ptrdiff_t movingBin = baseBin + offset; + if (movingBin < 0 || movingBin >= static_cast(bins)) { + continue; + } + const double weight = cubicBSpline( + movingPosition - static_cast(movingBin)); + if (weight <= 0.0) continue; + jointHistogram[ + fixedBin * bins + static_cast(movingBin) + ] += weight; + } +} + +double mutualInformationFromHistogram( + const std::vector& jointHistogram, + std::size_t bins) { + const double nan = std::numeric_limits::quiet_NaN(); + const double total = std::accumulate( + jointHistogram.begin(), jointHistogram.end(), 0.0); + if (!(total > 0.0) || !std::isfinite(total)) return nan; + + std::vector px(bins, 0.0); + std::vector py(bins, 0.0); + + for (std::size_t fixedBin = 0; fixedBin < bins; ++fixedBin) { + for (std::size_t movingBin = 0; movingBin < bins; ++movingBin) { + const double pxy = jointHistogram[ + fixedBin * bins + movingBin + ] / total; + px[fixedBin] += pxy; + py[movingBin] += pxy; + } + } + + double mi = 0.0; + for (std::size_t fixedBin = 0; fixedBin < bins; ++fixedBin) { + for (std::size_t movingBin = 0; movingBin < bins; ++movingBin) { + const double pxy = jointHistogram[ + fixedBin * bins + movingBin + ] / total; + const double pxPy = px[fixedBin] * py[movingBin]; + if (pxy > 0.0 && pxPy > 0.0) { + mi += pxy * std::log(pxy / pxPy); + } + } + } + return mi; +} + +double mattesMiFromValues( + const Pixel* fixedValues, + const Pixel* movingValues, + std::size_t count, + std::size_t bins, + std::optional fixedRange = std::nullopt, + std::optional movingRange = std::nullopt) { + const double nan = std::numeric_limits::quiet_NaN(); + if (count != 0U && (fixedValues == nullptr || movingValues == nullptr)) { + throw std::invalid_argument("Input value pointer is null."); + } + if (count < 2U) return nan; + if (bins < 4U) { + throw std::invalid_argument( + "bins must be >= 4 for cubic B-spline smoothing."); + } + if (bins > std::numeric_limits::max() / bins) { + throw std::length_error("Histogram dimensions are too large."); + } + + if (!fixedRange.has_value() || !movingRange.has_value()) { + std::array fixedCounts{}; + std::array movingCounts{}; + for (std::size_t i = 0; i < count; ++i) { + ++fixedCounts[fixedValues[i]]; + ++movingCounts[movingValues[i]]; + } + if (!fixedRange.has_value()) { + fixedRange = minMaxRangeFromCounts(fixedCounts, count); + } + if (!movingRange.has_value()) { + movingRange = minMaxRangeFromCounts(movingCounts, count); + } + } + + if (!isValidRange(*fixedRange) || !isValidRange(*movingRange)) return nan; + + std::vector jointHistogram(bins * bins, 0.0); + for (std::size_t i = 0; i < count; ++i) { + addMattesPair( + static_cast(fixedValues[i]), + static_cast(movingValues[i]), + *fixedRange, + *movingRange, + bins, + jointHistogram); + } + return mutualInformationFromHistogram(jointHistogram, bins); +} + +} // namespace + +cv::Mat1d MatteMIMap( + const cv::Mat& fixed, + const cv::Mat& moving, + const cv::Mat& mask, + int bins) { + validateInputs(fixed, moving, mask, bins); + + constexpr ChunkSize chunkSize{}; + constexpr std::size_t minValidPixels = 100U; + constexpr double minValidFraction = 0.10; + constexpr double lowerPercentile = 0.5; + constexpr double upperPercentile = 99.5; + + const GlobalCounts global = collectGlobalCounts(fixed, moving, mask); + if (global.validPairs == 0U) { + throw std::invalid_argument("The mask contains no valid pixels."); + } + + const IntensityRange fixedRange{ + percentileFromCounts(global.fixed, global.validPairs, lowerPercentile), + percentileFromCounts(global.fixed, global.validPairs, upperPercentile) + }; + const IntensityRange movingRange{ + percentileFromCounts(global.moving, global.validPairs, lowerPercentile), + percentileFromCounts(global.moving, global.validPairs, upperPercentile) + }; + + if (!isValidRange(fixedRange)) { + throw std::invalid_argument("Invalid fixed intensity range."); + } + if (!isValidRange(movingRange)) { + throw std::invalid_argument("Invalid moving intensity range."); + } + + const int nRows = ceilDividePositive(fixed.rows, chunkSize.height); + const int nCols = ceilDividePositive(fixed.cols, chunkSize.width); + cv::Mat1d nmiMap(nRows, nCols); + nmiMap.setTo(cv::Scalar(std::numeric_limits::quiet_NaN())); + + const std::size_t maxChunkPixels = + static_cast(std::min(chunkSize.height, fixed.rows)) * + static_cast(std::min(chunkSize.width, fixed.cols)); + + std::vector fixedChunkValues; + std::vector movingChunkValues; + fixedChunkValues.reserve(maxChunkPixels); + movingChunkValues.reserve(maxChunkPixels); + + for (int row = 0; row < nRows; ++row) { + for (int col = 0; col < nCols; ++col) { + const int y0 = row * chunkSize.height; + const int x0 = col * chunkSize.width; + const int y1 = y0 + std::min(chunkSize.height, fixed.rows - y0); + const int x1 = x0 + std::min(chunkSize.width, fixed.cols - x0); + const std::size_t totalPixels = + static_cast(y1 - y0) * + static_cast(x1 - x0); + + fixedChunkValues.clear(); + movingChunkValues.clear(); + + for (int y = y0; y < y1; ++y) { + const Pixel* fixedRow = fixed.ptr(y); + const Pixel* movingRow = moving.ptr(y); + const Pixel* maskRow = mask.empty() ? nullptr : mask.ptr(y); + for (int x = x0; x < x1; ++x) { + if (maskRow != nullptr && maskRow[x] == 0U) continue; + fixedChunkValues.push_back(fixedRow[x]); + movingChunkValues.push_back(movingRow[x]); + } + } + + const std::size_t validPixels = fixedChunkValues.size(); + if (validPixels < minValidPixels) continue; + + const double validFraction = + static_cast(validPixels) / + static_cast(totalPixels); + if (validFraction < minValidFraction) continue; + + nmiMap(row, col) = mattesMiFromValues( + fixedChunkValues.data(), + movingChunkValues.data(), + validPixels, + static_cast(bins), + std::optional{fixedRange}, + std::optional{movingRange}); + } + } + + return nmiMap; +} + +cv::Mat1d chunkedMatteMIMap( + const cv::Mat& fixed, + const cv::Mat& moving, + const cv::Mat& mask, + int bins) { + return MatteMIMap(fixed, moving, mask, bins); +} + +double MatteMI( + const cv::Mat& fixed, + const cv::Mat& moving, + const cv::Mat& mask, + int bins) { + validateInputs(fixed, moving, mask, bins); + + const std::size_t binCount = static_cast(bins); + if (binCount > std::numeric_limits::max() / binCount) { + throw std::length_error("Histogram dimensions are too large."); + } + + const GlobalCounts global = collectGlobalCounts(fixed, moving, mask); + if (global.validPairs < 2U) { + return std::numeric_limits::quiet_NaN(); + } + + const IntensityRange fixedRange = + minMaxRangeFromCounts(global.fixed, global.validPairs); + const IntensityRange movingRange = + minMaxRangeFromCounts(global.moving, global.validPairs); + if (!isValidRange(fixedRange) || !isValidRange(movingRange)) { + return std::numeric_limits::quiet_NaN(); + } + + std::vector jointHistogram(binCount * binCount, 0.0); + for (int y = 0; y < fixed.rows; ++y) { + const Pixel* fixedRow = fixed.ptr(y); + const Pixel* movingRow = moving.ptr(y); + const Pixel* maskRow = mask.empty() ? nullptr : mask.ptr(y); + for (int x = 0; x < fixed.cols; ++x) { + if (maskRow != nullptr && maskRow[x] == 0U) continue; + addMattesPair( + static_cast(fixedRow[x]), + static_cast(movingRow[x]), + fixedRange, + movingRange, + binCount, + jointHistogram); + } + } + + return mutualInformationFromHistogram(jointHistogram, binCount); +} diff --git a/src/matte_mi.h b/src/matte_mi.h new file mode 100644 index 00000000..9850d14a --- /dev/null +++ b/src/matte_mi.h @@ -0,0 +1,24 @@ +#ifndef VOLTRON_MATTE_MI_H +#define VOLTRON_MATTE_MI_H + +#include + +cv::Mat1d MatteMIMap( + const cv::Mat& fixed, + const cv::Mat& moving, + const cv::Mat& mask, + int bins = 50); + +cv::Mat1d chunkedMatteMIMap( + const cv::Mat& fixed, + const cv::Mat& moving, + const cv::Mat& mask, + int bins = 50); + +double MatteMI( + const cv::Mat& fixed, + const cv::Mat& moving, + const cv::Mat& mask, + int bins = 50); + +#endif \ No newline at end of file diff --git a/src/metrics.cpp b/src/metrics.cpp new file mode 100644 index 00000000..2117d261 --- /dev/null +++ b/src/metrics.cpp @@ -0,0 +1,372 @@ +#include + +// OpenCV +#include +#include "opencv2/features2d.hpp" +#include "opencv2/shape/shape_transformer.hpp" + +// Internal functions +#include "auxiliary.h" +#include "image.h" +#include "matte_mi.h" + +// Namespaces +using namespace Rcpp; +using namespace std; +using namespace cv; + +//// +// Quality Control +//// + +// check distribution of registered points +double checkMappedGridDistribution(Mat &im, Mat &h){ + + // message + std::string message; + + // get image shape + int height = im.rows; + int width = im.cols; + int height_interval = height > 50 ? (double) height/50.0 : 1; + int width_interval = width > 50 ? (double) width/50.0 : 1; + + // perspective transformation of grid points + std::vector gridpoints; + for (double i = 0.0; i <= height; i += height_interval) { + for (double j = 0.0; j <= width; j += width_interval) { + gridpoints.push_back(cv::Point2f(j,i)); + } + } + + // register grid points + std::vector gridpoints_reg; + if (h.rows == 2){ + cv::transform(gridpoints, gridpoints_reg, h); + } else if(h.rows == 3) { + cv::perspectiveTransform(gridpoints, gridpoints_reg, h); + } + + // Compute the standard deviation of the transformed points + return cppSD(gridpoints_reg); +} + +bool checkMaskAbundance(Mat &mask){ + int j=0; + for (int i = 0; i < mask.rows; i++) { + if (mask.at(i)) { + j++; + } + } + return j > 6; +} + +// compare the distance between two sets of match points +double medianMappingDistance(std::vector &keypoints1, std::vector &keypoints2, Mat &h) { + std::vector keypoints1_warped; + if(keypoints1.size() > 0){ + if (h.rows == 2){ + cv::transform(keypoints1, keypoints1_warped, h); + } else { + cv::perspectiveTransform(keypoints1, keypoints1_warped, h); + } + } + + return medianDistances(keypoints1_warped, keypoints2); +} + +// calculate inlier percentage +int checkInlierPercentage(Mat &mask){ + int j=0; + for (int i = 0; i < mask.rows; i++) { + if (mask.at(i)) { + j++; + } + } + double ratio = (double) j/mask.rows; + double perc = round(100.0 * ratio); + return (int) perc; +} + +void maskKeypoints(std::vector &keypoints1_good, std::vector &keypoints2_good, + std::vector &keypoints1_masked, std::vector &keypoints2_masked, + std::vector &top_matches, Mat &mask) +{ + int j=0; + for (int i = 0; i < mask.rows; i++) { + if (mask.at(i)) { + keypoints1_masked.push_back(keypoints1_good[i]); + keypoints2_masked.push_back(keypoints2_good[i]); + top_matches.push_back(cv::DMatch(static_cast(j), static_cast(j), 0)); + j++; + } + } +} + +// check if keypoints are degenerate +bool checkDegenerate(double pts1, double pts2) { + + // get warning message + bool is_degenerate = FALSE; + if(pts1 < 1.0 || pts2 < 1.0){ + is_degenerate = TRUE; + Rcout << "WARNING: points may be in a degenerate configuration." << endl; + } + + return is_degenerate; +} + +cv::Mat generateOverlapMask(cv::Size dsize, + cv::Mat& h, + cv::Size ssize) +{ + // generate mask + cv::Mat mask = cv::Mat::ones(ssize, CV_8UC1) * 255; + cv::Mat warped; + + // Keep masks crisp: nearest-neighbor only. + const int interp = cv::INTER_NEAREST; + const int borderMode = cv::BORDER_CONSTANT; + const cv::Scalar borderValue(0); + + // warp mask + if (h.rows == 2){ + cv::warpAffine(mask, warped, h, dsize, + interp, borderMode, borderValue); + } else { + cv::warpPerspective(mask, warped, h, dsize, + interp, borderMode, borderValue); + } + + // Force binary mask again. + cv::threshold(warped, warped, 0, 255, cv::THRESH_BINARY); + return warped; +} + +cv::Mat generateOverlapMask(cv::Mat& ref_image, + Ptr& tps, + cv::Size ssize) +{ + // generate mask + cv::Mat mask = cv::Mat::ones(ssize, CV_8UC1) * 255; + + // Keep masks crisp: nearest-neighbor only. + const int interp = cv::INTER_NEAREST; + // const int borderMode = cv::BORDER_CONSTANT; + // const cv::Scalar borderValue(0); + + mask = warpTPSImage(ref_image, mask, tps, + ref_image.rows, ref_image.cols, interp); + + // Force binary mask again. + cv::threshold(mask, mask, 0, 255, cv::THRESH_BINARY); + return mask; +} + +// [[Rcpp::export]] +Rcpp::IntegerVector generateOverlapMask(Rcpp::NumericVector& dsize, + Rcpp::NumericMatrix& trans_mat, + Rcpp::NumericVector& ssize){ + cv::Mat h = numericMatrixToMat(trans_mat); + cv::Mat mask = generateOverlapMask(cv::Size((int) dsize[0], (int) dsize[1]), + h, + cv::Size((int) ssize[0], (int) ssize[1])); + return matToMask(mask); + // return matToImage(mask); +} + +double Entropy(cv::Mat& im1, cv::Mat& overlapMask, int bins = 256) { + + // Histogram settings + int histSize = 256; + float range[] = {0.0, 256.0}; + const float* histRange = {range}; + int channels[] = {0}; + + // Compute histograms + cv::Mat hist; + cv::calcHist(&im1, 1, channels, overlapMask, + hist, 1, &histSize, &histRange); + + // Normalize histograms + cv::normalize(hist, hist, 0, 1, cv::NORM_MINMAX); + + // Convert counts to probabilities + hist /= cv::sum(hist)[0]; + + double entropy = 0.0; + for (int r = 0; r < hist.rows; ++r) + { + const float* ptr = hist.ptr(r); + + for (int c = 0; c < hist.cols; ++c) + { + double p = ptr[c]; + + if (p > 0.0) + entropy -= p * std::log(p); + } + } + + return entropy; +} + +double jointEntropy(cv::Mat& im1, cv::Mat& im2, + cv::Mat& overlapMask, int bins = 256) { + + // 2D histogram parameters + int histSize[] = {bins, bins}; + float range[] = {0.f, 256.f}; + const float* ranges[] = {range, range}; + int channels[] = {0, 1}; + + // calculate histogram + cv::Mat images[] = {im1, im2}; + cv::Mat hist; + cv::calcHist(images, + 2, + channels, + overlapMask, + hist, + 2, + histSize, + ranges, + true, + false); + cv::normalize(hist, hist, 0, 1, cv::NORM_MINMAX); + + // Convert counts to probabilities + hist /= cv::sum(hist)[0]; + + double entropy = 0.0; + for (int r = 0; r < hist.rows; ++r) + { + const float* ptr = hist.ptr(r); + + for (int c = 0; c < hist.cols; ++c) + { + double p = ptr[c]; + + if (p > 0.0) + entropy -= p * std::log(p); + } + } + + return entropy; +} + +double MutualInfo(cv::Mat& im1, cv::Mat& im2, + cv::Mat& overlapMask, int bins = 256) { + double ent1=Entropy(im1, overlapMask, bins); + double ent2=Entropy(im2, overlapMask, bins); + double ent12=jointEntropy(im1, im2, overlapMask, bins); + return ent1+ent2-ent12; +} + +double NormalizedMutualInfo(cv::Mat& im1, cv::Mat& im2, + cv::Mat& overlapMask, int bins = 256) { + double ent1=Entropy(im1, overlapMask, bins); + double ent2=Entropy(im2, overlapMask, bins); + double ent12=jointEntropy(im1, im2, overlapMask, bins); + return (ent1+ent2)/ent12; +} + +std::map getAlignmentMetrics(Mat &im1, Mat &im2, + Mat &mask, std::string type){ + + // Metrics + std::map metrics; + + // Compute histograms + int histSize = 256; + float range[] = {0.0, 256.0}; + const float* histRange = {range}; + int channels[] = {0}; + cv::Mat hist1, hist2; + cv::calcHist(&im1, 1, channels, mask, + hist1, 1, &histSize, &histRange); + cv::calcHist(&im2, 1, channels, mask, + hist2, 1, &histSize, &histRange); + + // Normalize histograms + // cv::normalize(hist1, hist1, 0, 1, cv::NORM_MINMAX); + // cv::normalize(hist2, hist2, 0, 1, cv::NORM_MINMAX); + hist1 /= cv::sum(hist1)[0]; + hist2 /= cv::sum(hist2)[0]; + + // Summary + Rcout << "Alignment Accuracy (" << type << "): " << endl; + metrics["Intersection"] = cv::compareHist(hist1, hist2, cv::HISTCMP_INTERSECT); + metrics["Bhattacharyya"] = cv::compareHist(hist1, hist2, cv::HISTCMP_BHATTACHARYYA); + metrics["Matte's MI"] = MatteMI(im2, im1, mask, 50); + + Rcout << " Intersection: " << metrics["Intersection"] << std::endl; + Rcout << " Bhattacharyya: " << metrics["Bhattacharyya"] << std::endl; + Rcout << " Matte's MI: " << metrics["Matte's MI"] << std::endl; + + // old metrics, keep for comparison + //metrics.push_back(cv::compareHist(hist1, hist2, cv::HISTCMP_CHISQR)); + // metrics.push_back(jointEntropy(im1, im2, mask, histSize)); + // metrics.push_back(MutualInfo(im1, im2, mask, histSize)); + // metrics.push_back(NormalizedMutualInfo(im1, im2, mask, histSize)); + + return metrics; +} + +// do overall checks on keypoints and images +std::map getKeypointMetrics(std::vector &points1, + std::vector &points2, + Mat &im1, Mat &im2, + Mat &h, Mat &mask) { + + // metrics list + std::map metrics; + + // Alignment report + Rcout << "Keypoint Report: " << endl; + + // Report final keypoints + Rcout << " Calculated transformation matrix with " << points1.size() << " keypoints" << endl; + metrics["#Keypoints"] = points1.size(); + + // get inlier percentages + double ratio = checkInlierPercentage(mask); + Rcout << " Inlier Percentage: " << ratio << endl; + metrics["Inlier Perc."] = ratio; + + // points stand. dev. + double points1_sd = cppSD(points1); + double points2_sd = cppSD(points2); + Rcout << " Std dev of points: x=" << points1_sd << " y=" << points2_sd << endl; + metrics["sd query kpts (>1?)"] = points1_sd; + metrics["sd ref. kpts (>1?)"] = points2_sd; + + // degenerate ? + bool degenerate_points = checkDegenerate(points1_sd, points2_sd); + metrics["Degenerate"] = (double) degenerate_points; + + // check distribution of points + double stddev = checkMappedGridDistribution(im1, h); + Rcout << " Std dev of registered points: " << stddev << endl; + if(stddev < 1.0 || stddev > max(im2.rows, im2.cols)){ + Rcout << " WARNING: Transformation may be poor - transformed points grid seem to be concentrated!" << endl; + metrics["Degenerate"] = 1.0; + } + metrics["sd grid (in [w,h]?)"] = stddev; + + // warp keypoints and check median distances + double md = medianMappingDistance(points1, points2, h); + Rcout << " Median distance between points: " << md << endl; + if(md > 3){ + Rcout << " WARNING: Transformation may be poor - mean euclidean distance of mapped source and destination key points is high!" << endl; + } + metrics["Median distance"] = md; + + // report degenerate + if((bool) metrics["Degenerate"]){ + Rcout << " WARNING: Registration is degenerate!" << endl; + } + + // return is_degenerate; + return metrics; +} \ No newline at end of file diff --git a/src/metrics.h b/src/metrics.h new file mode 100644 index 00000000..be75a2ac --- /dev/null +++ b/src/metrics.h @@ -0,0 +1,52 @@ +#include "Rcpp.h" +#include +#include "opencv2/shape/shape_transformer.hpp" + +// Namespaces +using namespace Rcpp; +using namespace std; +using namespace cv; + +#ifndef METRICS_H +#define METRICS_H + +// check distribution of registered points +double checkMappedGridDistribution(cv::Mat &im, cv::Mat &h); + +bool checkMaskAbundance(cv::Mat &mask); + +// compare the distance between two sets of match points +double medianMappingDistance(std::vector &keypoints1, std::vector &keypoints2, cv::Mat &h); + +// calculate inlier percentage +int checkInlierPercentage(cv::Mat &mask); + +void maskKeypoints(std::vector &keypoints1_good, std::vector &keypoints2_good, + std::vector &keypoints1_masked, std::vector &keypoints2_masked, + std::vector &top_matches, cv::Mat &mask); + +// check if keypoints are degenerate +bool checkDegenerate(double pts1, double pts2); + +// generate overlap mask for alignment +cv::Mat generateOverlapMask(cv::Size dsize, + cv::Mat& h, + cv::Size ssize); + +cv::Mat generateOverlapMask(cv::Mat& ref_image, + Ptr& tps, + cv::Size ssize); + +// get alignment metrics +std::map getAlignmentMetrics(cv::Mat &im1, + cv::Mat &im2, + cv::Mat &mask, + std::string type); + +// do overall checks on keypoints and metrics +std::map getKeypointMetrics(std::vector &points1, + std::vector &points2, + cv::Mat &im1, cv::Mat &im2, + cv::Mat &h, cv::Mat &mask); + +#endif \ No newline at end of file