diff --git a/DESCRIPTION b/DESCRIPTION index d51056f..19c26de 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,9 +2,10 @@ Package: velox Type: Package Title: Fast Raster Manipulation and Extraction Version: 0.2.0.9002 -Date: 2018-02-13 +Date: 2019-02-12 Author: Philipp Hunziker Maintainer: Philipp Hunziker +Contributor: Matthias Weigand BugReports: https://github.com/hunzikp/velox/issues Description: C++ accelerated raster manipulation and extraction. Depends: R (>= 3.0.2) @@ -19,8 +20,8 @@ Imports: sf (>= 0.5-5) LinkingTo: Rcpp, BH (>= 1.65.0-1) NeedsCompilation: yes -SystemRequirements: C++11 +SystemRequirements: C++11 Packaged: 2017-11-30 12:15 UTC; hunzikp LazyData: TRUE -RoxygenNote: 5.0.1 +RoxygenNote: 6.1.1 Suggests: testthat diff --git a/R/RcppExports.R b/R/RcppExports.R index 6e31020..6a357a4 100755 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -21,6 +21,10 @@ meanfocal_cpp <- function(rasterband, weights, wrow, wcol, band) { .Call('_velox_meanfocal_cpp', PACKAGE = 'velox', rasterband, weights, wrow, wcol, band) } +stddevfocal_cpp <- function(rasterband, wrow, wcol, band) { + .Call('_velox_stddevfocal_cpp', PACKAGE = 'velox', rasterband, wrow, wcol, band) +} + hittest_cpp <- function(rasterbands, dim, extent, res, polyX, polyY, polyCorners) { .Call('_velox_hittest_cpp', PACKAGE = 'velox', rasterbands, dim, extent, res, polyX, polyY, polyCorners) } diff --git a/R/velox_focal.R b/R/velox_focal.R index 17d5c02..b47f1a4 100755 --- a/R/velox_focal.R +++ b/R/velox_focal.R @@ -39,6 +39,47 @@ VeloxRaster$methods(medianFocal = function(wrow, wcol, bands=1) { +#' @title Standard Deviation focal +#' +#' @name VeloxRaster_stddevFocal +#' +#' @description +#' Applies a standard deviation filter of dimension \code{wcol x wrow} to a VeloxRaster. +#' +#' @details +#' Padding is currently not implemented. +#' +#' @param wrow y dimension of filter. Must be uneven integer. +#' @param wcol x dimension of filter. Must be uneven integer. +#' @param bands Numeric vector indicating bands where filter is applied. +#' +#' @return Void. +#' +#' @examples +#' ## Make VeloxRaster with two bands +#' mat1 <- matrix(1:100, 10, 10) +#' mat2 <- matrix(100:1, 10, 10) +#' vx <- velox(list(mat1, mat2), extent=c(0,1,0,1), res=c(0.1,0.1), +#' crs="+proj=longlat +datum=WGS84 +no_defs") +#' ## Median focal +#' vx$stddevFocal(wrow=5, wcol=5, bands=c(1,2)) +#' +NULL +VeloxRaster$methods(stddevFocal = function(wrow, wcol, bands=1) { + "See \\code{\\link{VeloxRaster_stddevFocal}}." + if (any(!(bands %in% 1:nbands))) { + stop(paste("VeloxRaster only has", nbands, "bands.")) + } + if (wrow < 0 | wcol < 0 | (wrow %% 2) == 0 | (wcol %% 2) == 0) { + stop(paste("wrow and wcol must be positive uneven integers.")) + } + for (i in bands) { + rasterbands[[i]] <<- stddevfocal_cpp(rasterband=rasterbands[[i]], wrow=wrow, wcol=wcol, band=i) + } +}) + + + #' @title Sum focal #' #' @name VeloxRaster_sumFocal diff --git a/man/BoostBoxGrid-class.Rd b/man/BoostBoxGrid-class.Rd index 29f276b..6e46372 100644 --- a/man/BoostBoxGrid-class.Rd +++ b/man/BoostBoxGrid-class.Rd @@ -7,4 +7,3 @@ \description{ This is a class for storing Rcpp pointers to C++ BoxGrid objects. } - diff --git a/man/BoostFactory.Rd b/man/BoostFactory.Rd index e3e4b5f..de6e596 100644 --- a/man/BoostFactory.Rd +++ b/man/BoostFactory.Rd @@ -6,4 +6,3 @@ \description{ Rcpp pointer to BoostFactory. } - diff --git a/man/BoostGeometries-class.Rd b/man/BoostGeometries-class.Rd index 83d6f4d..a72bdd6 100644 --- a/man/BoostGeometries-class.Rd +++ b/man/BoostGeometries-class.Rd @@ -7,4 +7,3 @@ \description{ This is a virtual class for storing Rcpp pointers to C++ GeometryCollection objects. } - diff --git a/man/BoostGrid-class.Rd b/man/BoostGrid-class.Rd index 5491f8c..5782cee 100644 --- a/man/BoostGrid-class.Rd +++ b/man/BoostGrid-class.Rd @@ -7,4 +7,3 @@ \description{ This is a virtual class for storing Rcpp pointers to C++ grid objects. } - diff --git a/man/BoostMultiLines-class.Rd b/man/BoostMultiLines-class.Rd index 22e3242..df4b900 100644 --- a/man/BoostMultiLines-class.Rd +++ b/man/BoostMultiLines-class.Rd @@ -7,4 +7,3 @@ \description{ This is a class for storing Rcpp pointers to C++ MultiLineCollection objects. } - diff --git a/man/BoostMultiPoints-class.Rd b/man/BoostMultiPoints-class.Rd index 877aaca..9662279 100644 --- a/man/BoostMultiPoints-class.Rd +++ b/man/BoostMultiPoints-class.Rd @@ -7,4 +7,3 @@ \description{ This is a class for storing Rcpp pointers to C++ MultiPointCollection objects. } - diff --git a/man/BoostMultiPolygons-class.Rd b/man/BoostMultiPolygons-class.Rd index 62b6cd1..a74d8d4 100644 --- a/man/BoostMultiPolygons-class.Rd +++ b/man/BoostMultiPolygons-class.Rd @@ -7,4 +7,3 @@ \description{ This is a class for storing Rcpp pointers to C++ MultiPolygonCollection objects. } - diff --git a/man/BoostPointGrid-class.Rd b/man/BoostPointGrid-class.Rd index 81d2901..b3b256a 100644 --- a/man/BoostPointGrid-class.Rd +++ b/man/BoostPointGrid-class.Rd @@ -7,4 +7,3 @@ \description{ This is a class for storing Rcpp pointers to C++ PointGrid objects. } - diff --git a/man/BoxGrid.Rd b/man/BoxGrid.Rd index 921255c..bda8f54 100644 --- a/man/BoxGrid.Rd +++ b/man/BoxGrid.Rd @@ -6,4 +6,3 @@ \description{ Rcpp pointer to BoxGrid. } - diff --git a/man/MultiLineCollection.Rd b/man/MultiLineCollection.Rd index b90f964..64a29df 100644 --- a/man/MultiLineCollection.Rd +++ b/man/MultiLineCollection.Rd @@ -6,4 +6,3 @@ \description{ Rcpp pointer to MultiLineCollection. } - diff --git a/man/MultiPointCollection.Rd b/man/MultiPointCollection.Rd index c511a50..f29fcf3 100644 --- a/man/MultiPointCollection.Rd +++ b/man/MultiPointCollection.Rd @@ -6,4 +6,3 @@ \description{ Rcpp pointer to MultiPointCollection. } - diff --git a/man/MultiPolygonCollection.Rd b/man/MultiPolygonCollection.Rd index edfb6bc..3afea92 100644 --- a/man/MultiPolygonCollection.Rd +++ b/man/MultiPolygonCollection.Rd @@ -6,4 +6,3 @@ \description{ Rcpp pointer to MultiPolygonCollection. } - diff --git a/man/PointGrid.Rd b/man/PointGrid.Rd index b00e2d5..ff1e1d6 100644 --- a/man/PointGrid.Rd +++ b/man/PointGrid.Rd @@ -6,4 +6,3 @@ \description{ Rcpp pointer to PointGrid. } - diff --git a/man/VeloxRaster-class.Rd b/man/VeloxRaster-class.Rd index 61ffc62..dedc303 100755 --- a/man/VeloxRaster-class.Rd +++ b/man/VeloxRaster-class.Rd @@ -2,8 +2,8 @@ % Please edit documentation in R/velox.R \docType{class} \name{VeloxRaster-class} -\alias{VeloxRaster} \alias{VeloxRaster-class} +\alias{VeloxRaster} \title{A Reference Class for velox rasters} \description{ A Reference Class for velox rasters @@ -23,6 +23,7 @@ A Reference Class for velox rasters \item{\code{crs}}{Coordinate reference system (Proj4 string)} }} + \section{Methods}{ \describe{ @@ -47,10 +48,10 @@ A Reference Class for velox rasters \item{\code{extract_points(sp)}}{See \code{\link{VeloxRaster_extract_points}}.} -\item{\code{getCoordinates()}}{See \code{\link{VeloxRaster_getCoordinates}}.} - \item{\code{get_data_type()}}{See \code{\link{VeloxRaster_get_data_type}}.} +\item{\code{getCoordinates()}}{See \code{\link{VeloxRaster_getCoordinates}}.} + \item{\code{im2col(wrow, wcol, band, padval = 0, rowframe = 0, colframe = 0, rowstride = 1, colstride = 1)}}{See \code{\link{VeloxRaster_im2col}}.} @@ -60,10 +61,13 @@ A Reference Class for velox rasters \item{\code{rasterize(spdf, field, band = 1, background = NULL, small = FALSE)}}{See \code{\link{VeloxRaster_rasterize}}.} +\item{\code{stddevFocal(wrow, wcol, bands = 1)}}{See \code{\link{VeloxRaster_stddevFocal}}.} + \item{\code{sumFocal(weights, bands = 1)}}{See \code{\link{VeloxRaster_sumFocal}}.} \item{\code{write(path, overwrite = FALSE)}}{See \code{\link{VeloxRaster_write}}.} }} + \examples{ ## Make VeloxRaster objects using the 'velox' function mat <- matrix(1:100, 10, 10) @@ -71,4 +75,3 @@ vx <- velox(mat, extent=c(0,1,0,1), res=c(0.1,0.1), crs="+proj=longlat +datum=WG class(vx) } - diff --git a/man/VeloxRaster_aggregate.Rd b/man/VeloxRaster_aggregate.Rd index 015e116..660457f 100755 --- a/man/VeloxRaster_aggregate.Rd +++ b/man/VeloxRaster_aggregate.Rd @@ -25,4 +25,3 @@ vx <- velox(mat, extent=c(0,1,0,1), res=c(0.1,0.1), crs="+proj=longlat +datum=WG ## Aggregate vx$aggregate(factor=c(2,2), aggtype='sum') } - diff --git a/man/VeloxRaster_as.RasterBrick.Rd b/man/VeloxRaster_as.RasterBrick.Rd index 9bb68f0..49ed180 100644 --- a/man/VeloxRaster_as.RasterBrick.Rd +++ b/man/VeloxRaster_as.RasterBrick.Rd @@ -24,4 +24,3 @@ library(raster) rs <- vx$as.RasterBrick() } - diff --git a/man/VeloxRaster_as.RasterLayer.Rd b/man/VeloxRaster_as.RasterLayer.Rd index 1d730ee..17a77cb 100755 --- a/man/VeloxRaster_as.RasterLayer.Rd +++ b/man/VeloxRaster_as.RasterLayer.Rd @@ -24,4 +24,3 @@ library(raster) rl <- vx$as.RasterLayer(band=1) } - diff --git a/man/VeloxRaster_as.RasterStack.Rd b/man/VeloxRaster_as.RasterStack.Rd index 6223b7e..df1a126 100755 --- a/man/VeloxRaster_as.RasterStack.Rd +++ b/man/VeloxRaster_as.RasterStack.Rd @@ -24,4 +24,3 @@ library(raster) rs <- vx$as.RasterStack() } - diff --git a/man/VeloxRaster_as.matrix.Rd b/man/VeloxRaster_as.matrix.Rd index 9ba5c93..d51d973 100755 --- a/man/VeloxRaster_as.matrix.Rd +++ b/man/VeloxRaster_as.matrix.Rd @@ -21,4 +21,3 @@ vx.mat <- vx$as.matrix(band=1) identical(mat, vx.mat) } - diff --git a/man/VeloxRaster_col2im.Rd b/man/VeloxRaster_col2im.Rd index ce60340..afe067c 100755 --- a/man/VeloxRaster_col2im.Rd +++ b/man/VeloxRaster_col2im.Rd @@ -44,4 +44,3 @@ vx$col2im(mat=patch.mat, wrow=2, wcol=2, band=1, rowframe=1, colframe=1, rowstri isTRUE(all.equal(mat, vx$as.matrix())) } - diff --git a/man/VeloxRaster_crop.Rd b/man/VeloxRaster_crop.Rd index fdf1d66..01e4d0f 100755 --- a/man/VeloxRaster_crop.Rd +++ b/man/VeloxRaster_crop.Rd @@ -24,4 +24,3 @@ vx$crop(c(0.3,0.7,0.3,0.7)) vx$extent } - diff --git a/man/VeloxRaster_drop.Rd b/man/VeloxRaster_drop.Rd index da40644..7a7e62a 100755 --- a/man/VeloxRaster_drop.Rd +++ b/man/VeloxRaster_drop.Rd @@ -22,4 +22,3 @@ vx <- velox(list(mat1, mat2), extent=c(0,1,0,1), res=c(0.1,0.1), vx$drop(bands=2) } - diff --git a/man/VeloxRaster_extract.Rd b/man/VeloxRaster_extract.Rd index 4a6d6ea..5ec4d00 100755 --- a/man/VeloxRaster_extract.Rd +++ b/man/VeloxRaster_extract.Rd @@ -49,4 +49,3 @@ spols <- gBuffer(spgeom=spoint, width=0.5) vx$extract(sp=spols, fun=mean) } - diff --git a/man/VeloxRaster_extract_points.Rd b/man/VeloxRaster_extract_points.Rd index 6a0d9f2..36ba044 100644 --- a/man/VeloxRaster_extract_points.Rd +++ b/man/VeloxRaster_extract_points.Rd @@ -28,4 +28,3 @@ spoint <- SpatialPoints(coords=coord) vx$extract_points(sp=spoint) } - diff --git a/man/VeloxRaster_getCoordinates.Rd b/man/VeloxRaster_getCoordinates.Rd index 6f43e1b..f56952b 100755 --- a/man/VeloxRaster_getCoordinates.Rd +++ b/man/VeloxRaster_getCoordinates.Rd @@ -17,4 +17,3 @@ vx <- velox(mat, extent=c(0,1,0,1), res=c(0.1,0.1), crs="+proj=longlat +datum=WG cmat <- vx$getCoordinates() } - diff --git a/man/VeloxRaster_get_data_type.Rd b/man/VeloxRaster_get_data_type.Rd index 5dcc441..d28f348 100644 --- a/man/VeloxRaster_get_data_type.Rd +++ b/man/VeloxRaster_get_data_type.Rd @@ -13,4 +13,3 @@ Get data type of a VeloxRaster Note that this method returns the data type of the raster, not the storage mode. Except in special cases, velox stores all raster data as double precision matrices. } - diff --git a/man/VeloxRaster_im2col.Rd b/man/VeloxRaster_im2col.Rd index 86e75fc..6fe3b0d 100755 --- a/man/VeloxRaster_im2col.Rd +++ b/man/VeloxRaster_im2col.Rd @@ -42,4 +42,3 @@ patch.mat <- vx$im2col(wrow=2, wcol=2, band=1, padval=0, dim(patch.mat) } - diff --git a/man/VeloxRaster_meanFocal.Rd b/man/VeloxRaster_meanFocal.Rd index 3debad4..e1e1a52 100755 --- a/man/VeloxRaster_meanFocal.Rd +++ b/man/VeloxRaster_meanFocal.Rd @@ -28,4 +28,3 @@ weights <- matrix(1, 5, 5) vx$meanFocal(weights=weights, bands=c(1,2)) } - diff --git a/man/VeloxRaster_medianFocal.Rd b/man/VeloxRaster_medianFocal.Rd index a43152c..b79479c 100755 --- a/man/VeloxRaster_medianFocal.Rd +++ b/man/VeloxRaster_medianFocal.Rd @@ -29,4 +29,3 @@ vx <- velox(list(mat1, mat2), extent=c(0,1,0,1), res=c(0.1,0.1), vx$medianFocal(wrow=5, wcol=5, bands=c(1,2)) } - diff --git a/man/VeloxRaster_rasterize.Rd b/man/VeloxRaster_rasterize.Rd index d462cf1..8cf12e2 100755 --- a/man/VeloxRaster_rasterize.Rd +++ b/man/VeloxRaster_rasterize.Rd @@ -45,4 +45,3 @@ spdf <- SpatialPolygonsDataFrame(Sr=spols, data=data.frame(value=1), match.ID=FA vx$rasterize(spdf=spdf, field="value", background=-1) } - diff --git a/man/VeloxRaster_stddevFocal.Rd b/man/VeloxRaster_stddevFocal.Rd new file mode 100644 index 0000000..1bb19ec --- /dev/null +++ b/man/VeloxRaster_stddevFocal.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/velox_focal.R +\name{VeloxRaster_stddevFocal} +\alias{VeloxRaster_stddevFocal} +\title{Standard Deviation focal} +\arguments{ +\item{wrow}{y dimension of filter. Must be uneven integer.} + +\item{wcol}{x dimension of filter. Must be uneven integer.} + +\item{bands}{Numeric vector indicating bands where filter is applied.} +} +\value{ +Void. +} +\description{ +Applies a standard deviation filter of dimension \code{wcol x wrow} to a VeloxRaster. +} +\details{ +Padding is currently not implemented. +} +\examples{ +## Make VeloxRaster with two bands +mat1 <- matrix(1:100, 10, 10) +mat2 <- matrix(100:1, 10, 10) +vx <- velox(list(mat1, mat2), extent=c(0,1,0,1), res=c(0.1,0.1), + crs="+proj=longlat +datum=WGS84 +no_defs") +## Median focal +vx$stddevFocal(wrow=5, wcol=5, bands=c(1,2)) + +} diff --git a/man/VeloxRaster_sumFocal.Rd b/man/VeloxRaster_sumFocal.Rd index 38dd178..413f71c 100755 --- a/man/VeloxRaster_sumFocal.Rd +++ b/man/VeloxRaster_sumFocal.Rd @@ -28,4 +28,3 @@ weights <- matrix(1, 5, 5) vx$sumFocal(weights=weights, bands=c(1,2)) } - diff --git a/man/VeloxRaster_write.Rd b/man/VeloxRaster_write.Rd index 73caae1..cb92758 100755 --- a/man/VeloxRaster_write.Rd +++ b/man/VeloxRaster_write.Rd @@ -14,4 +14,3 @@ Void. \description{ Write a VeloxRaster to disk as a GeoTiff file } - diff --git a/man/bg_intersects.generic.Rd b/man/bg_intersects.generic.Rd index 50d15b4..6c75f23 100644 --- a/man/bg_intersects.generic.Rd +++ b/man/bg_intersects.generic.Rd @@ -3,12 +3,12 @@ \docType{methods} \name{bg_intersects} \alias{bg_intersects} +\alias{bg_intersects,BoostMultiPolygons,BoostBoxGrid-method} \alias{bg_intersects,BoostMultiLines,BoostBoxGrid-method} +\alias{bg_intersects,BoostMultiPolygons,BoostPointGrid-method} +\alias{bg_intersects,BoostMultiPolygons,BoostGeometries-method} \alias{bg_intersects,BoostMultiLines,BoostGeometries-method} \alias{bg_intersects,BoostMultiPoints,BoostGeometries-method} -\alias{bg_intersects,BoostMultiPolygons,BoostBoxGrid-method} -\alias{bg_intersects,BoostMultiPolygons,BoostGeometries-method} -\alias{bg_intersects,BoostMultiPolygons,BoostPointGrid-method} \title{Test whether two BoostObjects intersect} \usage{ bg_intersects(obj1, obj2) @@ -45,4 +45,3 @@ list(rbind(c(0,0), c(2,0), c(2,2), c(0,2), c(0,0)))))) bg_intersects(pol, pts) } - diff --git a/man/boost.Rd b/man/boost.Rd index ccd1b1a..3897965 100644 --- a/man/boost.Rd +++ b/man/boost.Rd @@ -2,13 +2,13 @@ % Please edit documentation in R/boostgeom.R \name{boost} \alias{boost} -\alias{boost.VeloxRaster} -\alias{boost.sfc_LINESTRING} +\alias{boost.sfc_MULTIPOLYGON} +\alias{boost.sfc_POLYGON} \alias{boost.sfc_MULTILINESTRING} +\alias{boost.sfc_LINESTRING} \alias{boost.sfc_MULTIPOINT} -\alias{boost.sfc_MULTIPOLYGON} \alias{boost.sfc_POINT} -\alias{boost.sfc_POLYGON} +\alias{boost.VeloxRaster} \title{Cast a sfc object as a BoostObject} \usage{ boost(x, ...) @@ -47,4 +47,3 @@ sfc <- sf::st_sfc(sf::st_point(c(0,1))) boostPoints <- boost(sfc) } - diff --git a/man/length.BoostGeometries.Rd b/man/length.BoostGeometries.Rd index d027763..358845b 100644 --- a/man/length.BoostGeometries.Rd +++ b/man/length.BoostGeometries.Rd @@ -16,4 +16,3 @@ An integer scalar. \description{ Returns the length (number of Geometries) of a BoostGeometries object. } - diff --git a/man/plot.BoostGeometries.Rd b/man/plot.BoostGeometries.Rd index 77da601..08fc008 100644 --- a/man/plot.BoostGeometries.Rd +++ b/man/plot.BoostGeometries.Rd @@ -20,4 +20,3 @@ Void. \description{ Plot a BoostGeometries object using the \code{sf} plotting function. } - diff --git a/man/subset.BoostGeometries.Rd b/man/subset.BoostGeometries.Rd index 61230aa..10e8632 100644 --- a/man/subset.BoostGeometries.Rd +++ b/man/subset.BoostGeometries.Rd @@ -2,9 +2,9 @@ % Please edit documentation in R/boostgeom.R \docType{methods} \name{[,BoostMultiPolygons,ANY,ANY-method} +\alias{[,BoostMultiPolygons,ANY,ANY-method} \alias{[,BoostMultiLines,ANY,ANY-method} \alias{[,BoostMultiPoints,ANY,ANY-method} -\alias{[,BoostMultiPolygons,ANY,ANY-method} \title{Subset a BoostGeometries object} \usage{ \S4method{[}{BoostMultiPolygons,ANY,ANY}(x, i) @@ -24,4 +24,3 @@ A BoostGeometries object. \description{ Extract a subset of geometries from a BoostGeometries object. } - diff --git a/man/unboost.Rd b/man/unboost.Rd index e807348..d46069b 100644 --- a/man/unboost.Rd +++ b/man/unboost.Rd @@ -2,9 +2,9 @@ % Please edit documentation in R/boostgeom.R \name{unboost} \alias{unboost} +\alias{unboost.BoostMultiPolygons} \alias{unboost.BoostMultiLines} \alias{unboost.BoostMultiPoints} -\alias{unboost.BoostMultiPolygons} \title{Cast a BoostGeometries object as a sfc object} \usage{ unboost(x, ...) @@ -37,4 +37,3 @@ sfc2 <- unboost(boostPoints) print(identical(sfc, sfc2)) } - diff --git a/man/velox.Rd b/man/velox.Rd index cf14489..a0912ef 100755 --- a/man/velox.Rd +++ b/man/velox.Rd @@ -40,4 +40,3 @@ mat.ls <- list(mat1, mat2) vx <- velox(mat.ls, extent=c(0,1,0,1), res=c(0.1,0.1), crs="+proj=longlat +datum=WGS84 +no_defs") } - diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 3aac839..75c5dec 100755 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -77,6 +77,20 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// stddevfocal_cpp +NumericMatrix stddevfocal_cpp(NumericMatrix rasterband, int wrow, int wcol, int band); +RcppExport SEXP _velox_stddevfocal_cpp(SEXP rasterbandSEXP, SEXP wrowSEXP, SEXP wcolSEXP, SEXP bandSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< NumericMatrix >::type rasterband(rasterbandSEXP); + Rcpp::traits::input_parameter< int >::type wrow(wrowSEXP); + Rcpp::traits::input_parameter< int >::type wcol(wcolSEXP); + Rcpp::traits::input_parameter< int >::type band(bandSEXP); + rcpp_result_gen = Rcpp::wrap(stddevfocal_cpp(rasterband, wrow, wcol, band)); + return rcpp_result_gen; +END_RCPP +} // hittest_cpp NumericMatrix hittest_cpp(List rasterbands, NumericVector dim, NumericVector extent, NumericVector res, NumericVector polyX, NumericVector polyY, int polyCorners); RcppExport SEXP _velox_hittest_cpp(SEXP rasterbandsSEXP, SEXP dimSEXP, SEXP extentSEXP, SEXP resSEXP, SEXP polyXSEXP, SEXP polyYSEXP, SEXP polyCornersSEXP) { @@ -197,6 +211,7 @@ static const R_CallMethodDef CallEntries[] = { {"_velox_medianfocal_cpp", (DL_FUNC) &_velox_medianfocal_cpp, 4}, {"_velox_sumfocal_cpp", (DL_FUNC) &_velox_sumfocal_cpp, 5}, {"_velox_meanfocal_cpp", (DL_FUNC) &_velox_meanfocal_cpp, 5}, + {"_velox_stddevfocal_cpp", (DL_FUNC) &_velox_stddevfocal_cpp, 4}, {"_velox_hittest_cpp", (DL_FUNC) &_velox_hittest_cpp, 7}, {"_velox_unhit_cpp", (DL_FUNC) &_velox_unhit_cpp, 4}, {"_velox_im2col_cpp", (DL_FUNC) &_velox_im2col_cpp, 10}, diff --git a/src/focal.cpp b/src/focal.cpp index 753bd04..059efc4 100755 --- a/src/focal.cpp +++ b/src/focal.cpp @@ -1,6 +1,7 @@ #include #include #include "median.h" +#include "stddev.h" using namespace Rcpp; using namespace std; @@ -108,4 +109,34 @@ NumericMatrix meanfocal_cpp(NumericMatrix rasterband, NumericMatrix weights, int return(newband); } - +// [[Rcpp::export]] +NumericMatrix stddevfocal_cpp(NumericMatrix rasterband, int wrow, int wcol, int band) { + + int nrow = rasterband.nrow(); + int ncol = rasterband.ncol(); + + NumericMatrix newband(nrow, ncol); + vector v; + + int idim = (wrow-1)/2; + int jdim = (wcol-1)/2; + + for (int i = 0; i < nrow; i++) { + int imin = i - idim; + int imax = i + idim + 1; + for (int j = 0; j < ncol; j++) { + int jmin = j - jdim; + int jmax = j + jdim + 1; + for (int wi = imin; wi < imax; wi++) { + for (int wj = jmin; wj < jmax; wj++) { + if (wi >= 0 && wi < nrow && wj >= 0 && wj < ncol) { + v.push_back(rasterband(wi, wj)); + } + } + } + newband(i,j) = stddev(v); + v.clear(); + } + } + return(newband); +} diff --git a/src/integral.cpp b/src/integral.cpp new file mode 100644 index 0000000..18857b8 --- /dev/null +++ b/src/integral.cpp @@ -0,0 +1,237 @@ +#include +#include +#include +using namespace Rcpp; + +// [[Rcpp::export]] +NumericMatrix integralRaster(NumericMatrix rasterband) { + + int rows = rasterband.nrow(); + int cols = rasterband.ncol(); + + for (int r = 0; r < rows; r++){ + for (int c = 0; c < cols; c++){ + if (c == 0){ + if (r != 0){ + rasterband(r, c) = rasterband(r, c) + rasterband(r - 1, c); + } + } else { + if(r == 0){ + rasterband(r, c) = rasterband(r, c) + rasterband(r, c - 1); + } else { + rasterband(r, c) = rasterband(r, c) + rasterband(r - 1, c) + + rasterband(r, c - 1) - rasterband(r - 1, c - 1); + } + } + } + } + return rasterband; +} + + +void fillKernel(NumericMatrix kernel, int diameter){ + // ancilliary function which uses the Pythagorean theorem to determine which + // cell of the disc shaped focal kernel is within the disc shaped + + int radius = (diameter-1) / 2; + int centX = radius; + int centY = radius; + int dX, dY; + float dist; + + for (int y = 0; y < diameter; y++) { + dY = abs(y - centY); + for (int x = 0; x < diameter; x++){ + dX = abs(x - centX); + dist = sqrt(pow(dX, 2) + pow(dY, 2)); + if(dist <= radius) { + kernel(x,y) = 1; + } + } + } +} + +void fillFirstRows(int* firstRows, NumericMatrix kernel, int diameter){ + // group columns of kernel into boxes with same first and last row to make + // even more use of the integral image by creating rectangular boxes + + int prevFirstRow; + int currValue; + int prevValue; + + // iterate over cols to find first row of the kernel in each col; + for (int c = 0; c < diameter; c++){ + + // reset last value at beginning of each column + prevValue = 0; + + // go through rows and check whether the value changed + for (int r = 0; r < diameter; r++){ + + // get the current value in the kernel matrix + currValue = kernel(r,c); + + // if the value changed, then disc shape started/ended for this column + if (prevValue != currValue){ + + // if value changed from 0 to 1 + if (prevValue == 0){ + firstRows[c] = r; + if (r != prevFirstRow) { + prevFirstRow = r; + } + prevValue = currValue; + } + } + } + // return firstRows; + } +} + + + + +std::vector> getBoxCoordinates(int * firstRows, int nRows){ + // anciliary function to get the upper left and lower right coordinates of + // boxes. uses the first rows of the squared kernel which are occupied by the + // circular shape if the filter to create groups of rectangular boxes which + // start at the same line in the matrix. That is, the algorithm can make use + // of the integral image and gain in efficiency. + // This function returns a list of the upper left / lower right coordinates + // of the individual boxes + std::vector> ullr; + + int prevVal; // previous value of firstRows, needed for detecting change + int groupStartX = 0; // remember the column in which the current group started + std::vector box; // initialize vector for ullr coordinates of group + + + for (int x = 0; x <= nRows; x++) { + // at the first iteration set up group variable and remember the start of the + // group in x direction + if (x == 0){ + prevVal = firstRows[x]; + groupStartX = x; + } else { + // if the firstrow does not change (i.e. the kernel starts at the same + // line as in the column before), just go on. + + if (firstRows[x] != prevVal){ + // if the kernel changes the row, create a new box for the group (i.e. + // rectangle box) which has just finished. + + box.push_back(groupStartX); // left + box.push_back(prevVal); // upper + box.push_back(x - 1); // right + box.push_back(nRows - 1 - prevVal); // lower + + // add the box to the the output + ullr.push_back(box); + box.clear(); + + // reassign prevVal to currenct group's firstRow's + prevVal = firstRows[x]; + // and save x position + groupStartX = x; + } + } + } + + return ullr; +} + + +// [[Rcpp::export]] +NumericMatrix intRadialFocalMean(NumericMatrix rasterband, int diameter) { + // Implementation of a radial focal mean algorithm using itegral raster images. + // This implementaiton is advantagous to classical focal operators, as it is + // more efficient especially with large kernels. + // With radial filters, one probably has to sacrifice some efficiency, though. + + // create output raster band + int rows = rasterband.nrow(); + int cols = rasterband.ncol(); + NumericMatrix newband(rows, cols); + + // derive radius from diameter + int radius; + if (diameter % 2 == 1){ + radius = (diameter - 1) / 2; + } else { + radius = diameter / 2; + diameter = diameter + 1; + } + + // create disc shaped kernel + NumericMatrix kernel(diameter, diameter); + fillKernel(kernel, diameter); + + // count number of pixels in kernel; needed for creating the mean + int nPix = 0; + for (int i = 0; i < diameter; i++) + for (int j = 0; j < diameter; j++) + if (kernel(i,j) == 1) nPix++; + + // find first rows in each column which are part of the disc kernel + int firstRows[diameter]; + fillFirstRows(firstRows, kernel, diameter); + + // std::cout << "Firstrows: "; + // for (auto x : firstRows) std::cout << x << " "; + // std::cout << std::endl; + + // find rectangular boxes in the disc shaped kernel to speed up processing + std::vector> ullr = getBoxCoordinates(firstRows, diameter); + + // convert ullr coordinates to values relative to the centering pixel of the kernel + for (int x = 0; x < ullr.size(); x++) { + for (int y = 0; y < ullr[x].size(); y++) { + ullr[x][y] = ullr[x][y] - ((diameter - 1) / 2); + } + } + + + // from here on copied and adapted from original velox/src/focal.cpp + std::vector v; + int sum; + + int idim = (diameter-1)/2; + int jdim = (diameter-1)/2; + + for (int i = 0; i < rows; i++) { + int imin = i - idim; + int imax = i + idim + 1; + for (int j = 0; j < cols; j++) { + int jmin = j - jdim; + int jmax = j + jdim + 1; + + // for (int wi = imin; wi < imax; wi++) { + // for (int wj = jmin; wj < jmax; wj++) { + // if (wi >= 0 && wi < rows && wj >= 0 && wj < cols) { + // // v.push_back(rasterband(wi, wj)); + // } + // } + // } + + sum = accumulate(v.begin(),v.end(),0); + + + newband(i,j) = sum / nPix;; + v.clear(); + } + } + return kernel; + // return newband; +} + + +/*** R +m <- matrix(c(0), nrow = 9, ncol = 9, byrow = TRUE) +m[5,5] <- 1 + +(i <- integralRaster(m)) +(imn <- intRadialFocalMean(i, 9)) +raster::plot(raster::raster(imn)) +# largeM <- matrix(runif(100000000), nrow = 10000) +# largeMint <- integralRaster(largeM) +*/ diff --git a/src/stddev.cpp b/src/stddev.cpp new file mode 100644 index 0000000..a6d15e6 --- /dev/null +++ b/src/stddev.cpp @@ -0,0 +1,27 @@ +#include +#include + +using namespace std; + +double stddev(vector scores) { + + double stddev; + double sum = 0.0, mean, var = 0.0; + + size_t size = scores.size(); + + + for (int i = 0; i < size; i++){ + sum += scores[i]; + } + + mean = sum / size; + + for (int i = 0; i < size; i++){ + var += pow(scores[i] - mean, 2); + } + + stddev = sqrt(var / size); + + return stddev; +} diff --git a/src/stddev.h b/src/stddev.h new file mode 100644 index 0000000..fe04c3c --- /dev/null +++ b/src/stddev.h @@ -0,0 +1,7 @@ +#ifndef STDDEV_H_INCLUDED +#define STDDEV_H_INCLUDED + +#include +double stddev(std::vector scores); + +#endif