From f26bb26045ab7da0365f6fd1edf7339a8293f9f0 Mon Sep 17 00:00:00 2001
From: rhijmans <r.hijmans@gmail.com>
Date: Sat, 7 Oct 2023 01:22:51 -0700
Subject: [PATCH] fixes #1305

---
 R/generics.R   | 20 ++++++++++++++------
 man/stretch.Rd |  4 ++--
 2 files changed, 16 insertions(+), 8 deletions(-)

diff --git a/R/generics.R b/R/generics.R
index 61f501532..4ec9c94df 100644
--- a/R/generics.R
+++ b/R/generics.R
@@ -971,15 +971,23 @@ setMethod("scale", signature(x="SpatRaster"),
 
 
 setMethod("stretch", signature(x="SpatRaster"),
-	function(x, minv=0, maxv=255, minq=0, maxq=1, smin=NA, smax=NA, histeq=FALSE, scale=1, filename="", ...) {
-		opt <- spatOptions(filename, ...)
+	function(x, minv=0, maxv=255, minq=0, maxq=1, smin=NA, smax=NA, histeq=FALSE, scale=1, maxcell=500000, filename="", ...) {
 		if (histeq) {
-			eqStretch <- function(x){
-				ecdfun <- stats::ecdf(x)(x)
-				ecdfun(x)
+			if (nlyr(x) > 1) {
+				warn("stretch", "only the first layer of x is used")
+				x <- x[[1]]
 			}
-			setValues(x, apply(values(x), 2, function(i) stats::ecdf(i)(i))) * scale
+			if (scale == 1) {
+				ecdfun <- stats::ecdf(na.omit(spatSample(x, maxcell, "regular")[,1]))
+			} else {
+				ecdfun <- function(y) {
+					f <- stats::ecdf(na.omit(spatSample(x, maxcell, "regular")[,1]))
+					f(y) * scale
+				}
+			}
+			app(x, ecdfun, filename=filename, ...)
 		} else {
+			opt <- spatOptions(filename, ...)
 			x@cpp <- x@cpp$stretch(minv, maxv, minq, maxq, smin, smax, opt)
 			messages(x, "stretch")
 		}
diff --git a/man/stretch.Rd b/man/stretch.Rd
index 944a6b8cc..1bddbb25c 100644
--- a/man/stretch.Rd
+++ b/man/stretch.Rd
@@ -10,12 +10,12 @@ Linear or histogram equalization stretch of values in a SpatRaster.
 
 For linear stretch, provide the desired output range (\code{minv} and \code{maxv}) and the lower and upper bounds in the original data, either as quantiles (\code{minq} and \code{maxq}, or as cell values (\code{smin} and \code{smax}). If \code{smin} and \code{smax} are both not \code{NA}, \code{minq} and \code{maxq} are ignored.
 
-For histogram equalization, these arguments are ignored, but you can provide the desired scale of the output. 
+For histogram equalization, these arguments are ignored, but you can provide the desired scale of the output and the maximum number of cells that is used to compute the histogram (empirical cumulative distribution function). Only the first layer of the input SpatRaster is processed when using histogram equalization.
 }
 
 \usage{
 \S4method{stretch}{SpatRaster}(x, minv=0, maxv=255, minq=0, maxq=1, smin=NA, smax=NA,
-		histeq=FALSE, scale=1, filename="", ...)
+		histeq=FALSE, scale=1, maxcell=500000, filename="", ...)
 }
 
 \arguments{