Package 'rcaiman'

Description

Applies a method to each set of pixels defined by a direction and a constant field of view (FOV). By default, several built-in methods are available (see method), but a custom function can also be provided via the fun argument.

Usage

apply_by_direction(
  r,
  z,
  a,
  m,
  spacing = 10,
  laxity = 2.5,
  fov = c(30, 40, 50),
  method = c("thr_isodata", "thr_twocorner", "detect_bg_dn", "fit_coneshaped_model",
    "fit_trend_surface_np1", "fit_trend_surface_np6"),
  fun = NULL,
  parallel = TRUE,
  cores = NULL,
  logical = TRUE,
  leave_free = 1
)

Arguments

Value

terra::SpatRaster object with two layers: "dn" for digital number values and "n" for the number of valid pixels used in each directional estimate.

Note

This function is part of a manuscript currently under preparation.

References

There are no references for Rd macro ⁠\insertAllCites⁠ on this help page.

Examples

## Not run: 
caim <- read_caim()
r <- caim$Blue
z <- zenith_image(ncol(caim), lens())
a <- azimuth_image(z)
m <- !is.na(z)

# Automatic sky brightness estimation
sky <- apply_by_direction(r, z, a, m, spacing = 10, fov = c(30, 60),
                          method = "detect_bg_dn", parallel = TRUE)
plot(sky$dn)
plot(r / sky$dn)

# Using cone-shaped model
sky_cs <- apply_by_direction(caim, z, a, m, spacing = 15, fov = 60,
                             method = "fit_coneshaped_model", parallel = TRUE)
plot(sky_cs$dn)

# Using trend surface model
sky_s <- apply_by_direction(caim, z, a, m, spacing = 15, fov = 60,
                            method = "fit_trend_surface_np1", parallel = TRUE)
plot(sky_s$dn)

# Using a custom thresholding function
thr <- apply_by_direction(r, z, a, m, 15, fov = c(30, 40, 50),
  fun = function(r, z, a, m) {
    thr <- tryCatch(thr_twocorner(r[m])$tm, error = function(e) NA)
    r[] <- thr
    r
  },
  parallel = TRUE
)
plot(thr$dn)
plot(binarize_with_thr(r, thr$dn))

## End(Not run)

Description

Evaluate a user-supplied function locally using neighboring points defined on the hemispherical domain.

Usage

apply_locally(
  sampling_points,
  query_points = NULL,
  z,
  a,
  k = 20,
  angular_radius = 20,
  rule = "any",
  fun,
  parallel = TRUE,
  cores = NULL,
  logical = TRUE,
  leave_free = 1
)

Arguments

Details

The query_points argument lets the user evaluate fun at a custom set of locations. For each point in query_points, neighbors are searched within sampling_points. When query_points is NULL (default), the evaluation is performed for every point in sampling_points. Neighbor selection is based on spherical distance, see calc_spherical_distance(). The function fun is then applied to the corresponding subset of neighboring sampling points.

The input passed to fun is a subset of rows of sampling_points, therefore, fun must extract and process the column(s) of interest, e.g.: function(df) mean(df$dn).

The angular neighborhood is defined by two parameters:

• k: number of nearest neighbors used by fun.

• angular_radius: maximum allowed spherical distance (deg).

Argument rule controls how neighbor availability is interpreted:

• "any": all neighbors found within angular_radius are passed to fun (even if fewer than k).

• "all": if fewer than k neighbors lie within the radius, the query point is assigned NA.

Value

A data.frame equal to query_points with one additional column, named "out", containing the result of applying fun at each query point.

See Also

skygrid_centers() and fibonacci_points() for producing query_points.

Examples

## Not run: 
caim <- read_caim()
r <- caim$Blue
z <- zenith_image(ncol(caim), lens())
a <- azimuth_image(z)

# local thresholding
sampling_points <- fibonacci_points(z, a, 3)
display_caim(caim, sampling_points = sampling_points)
sampling_points <- extract_dn(r, sampling_points, use_window = FALSE)
head(sampling_points)
eval_points <- apply_locally(
  sampling_points,
  fibonacci_points(z, a, 10),
  z,
  a,
  k = 500,
  angular_radius = 30,
  rule = "any",
  fun = function(df) thr_isodata(df$Blue)
)

thr <- triangulate(eval_points[!is.na(eval_points$out),], r, col_id = 3)
plot(thr)
plot(binarize_with_thr(r, thr))

# Detect background DN
sampling_points <- fibonacci_points(z, a, 2)
sampling_points <- extract_dn(r, sampling_points, use_window = FALSE)
query_points <- fibonacci_points(z, a, 10)
eval_points <- apply_locally(
  sampling_points,
  query_points,
  z,
  a,
  k = 1000,
  angular_radius = 40,
  rule = "any",
  fun = function(df) {
    tryCatch(thr_twocorner(df[,3])$up, error = function(e) NA)
  },
  parallel = TRUE
)
i <- !is.na(eval_points$out)

sky <- triangulate(eval_points[i, ], r, col_id = 3)
plot(sky)
plot(r/sky)


## End(Not run)

Description

Quantify how a set of sky sampling points covers an equal-area hemispherical segmentation using (i) Kullback–Leibler divergence from a uniform distribution and (ii) coverage, defined as the fraction of cells receiving at least one point.

Usage

assess_sampling_uniformity(sky_points, equalarea_seg)

Arguments

Details

Points are projected onto an equal-area segmentation (equalarea_seg), and cell-wise counts are computed.

Let be $N_{\text{points}}$ the total number of points, $N_{\text{cells}}$ the total number of cells, and $c_i$ the point count in cell $i$ ($i = 1,\dots,N_{\text{cells}}$). The empirical distribution for cell $i$ is:

$p_i = \frac{c_i}{N_{\text{points}}}, \quad \text{with} \quad \sum_{i=1}^{N_{\text{cells}}} p_i = 1.$

The uniform reference is

$u_i = 1 / N_{\text{cells}}.$

Kullback–Leibler divergence is computed excluding zero-probability terms:

$\mathrm{KL}(p\|u) = \sum_{i: p_i > 0}^{N_{\text{cells}}} p_i \cdot \log(p_i \cdot N_{\text{cells}}).$

Coverage is defined as:

$\text{coverage} = \frac{\text{Count}(\{i : c_i > 0\})}{N_{\text{cells}}}.$

Value

List with two components:

kl

numeric value. Kullback–Leibler divergence from uniformity.

coverage

numeric value in $[0,1]$. Fraction of cells with at least one point.

Examples

## Not run: 
caim <- read_caim()
r <- caim$Blue
z <- zenith_image(ncol(caim), lens())
a <- azimuth_image(z)
m <- !is.na(z)

seg <- equalarea_segmentation(z, a, 20)
bin <- binarize_by_region(r, seg, method = "thr_twocorner")
bin <- bin & select_sky_region(z, 0, 80)

g <- skygrid_segmentation(z, a, 5)
sky_points <- sample_sky_points(r, bin, g, dist_to_black = 3)
plot(r)
points(sky_points$col, nrow(caim) - sky_points$row, col = 2, pch = 10)

seg <- equalarea_segmentation(z, a, 200)
kl <- lapply(1:10, function(i) {
  sky_points <-  rem_nearby_points(sky_points, NULL, z, a, i, space = "spherical")
  assess_sampling_uniformity(sky_points, seg)$kl
})
plot(unlist(kl), type = "l")

## End(Not run)

Description

Creates a single-layer raster in which pixel values represent azimuth angles, assuming an upwards-looking hemispherical photograph with the optical axis vertically aligned.

Usage

azimuth_image(z, orientation = 0)

Arguments

Value

terra::SpatRaster with the same dimensions as the input zenith image. Each pixel contains the azimuth angle in degrees, with zero representing North and angles increasing counter-clockwise. The object carries attributes orientation and lens_coef.

Note

If orientation = 0, North (0 deg) is located at the top of the image, as in conventional maps, but East (90 deg) and West (270 deg) appear flipped relative to maps. To understand this, take two flash-card-sized pieces of paper. Place one on a table in front of you and draw a compass rose on it. Hold the other above your head, with the side facing down toward you, and draw another compass rose following the directions from the one on the table. This mimics the situation of taking an upwards-looking photo with a smartphone while viewing the screen, and it will result in a mirrored arrangement. Compare both drawings to see the inversion.

Examples

z <- zenith_image(600, lens("Nikon_FCE9"))
a <- azimuth_image(z)
plot(a)
## Not run: 
a <- azimuth_image(z, 45)
plot(a)

## End(Not run)

Description

Perform thresholding of greyscale images by applying a method regionally, using a segmentation map.

Usage

binarize_by_region(r, seg, method)

Arguments

Details

This function supports several thresholding methods applied within the regions defined by seg:

Methods from the autothresholdr package:

Any method supported by autothresholdr::auto_thresh() can be used by specifying its name. For example, "IsoData" applies the classic iterative intermeans algorithm Ridler and Calvard (1978), which is among the most recommended for canopy photography (Jonckheere et al. 2005).

In-package implementation of IsoData:

Use "thr_isodata" to apply thr_isodata(), a native implementation of the same algorithm

Two-corner method:

Use "thr_twocorner" to apply thr_twocorner(), which implements a geometric thresholding strategy based on identifying inflection points in the histogram, first introduced to canopy photography by Macfarlane (2011). Since this method might fail, the fallback is thr_isodata

Value

Logical terra::SpatRaster (TRUE for sky, FALSE for non-sky) of the same dimensions as r.

Note

When methods from the autothresholdr package are used, r values should be constrained to the range $[0, 1]$. See normalize_minmax().

References

Jonckheere I, Nackaerts K, Muys B, Coppin P (2005). “Assessment of automatic gap fraction estimation of forests from digital hemispherical photography.” Agricultural and Forest Meteorology, 132(1-2), 96–114. doi:10.1016/j.agrformet.2005.06.003.

Leblanc SG, Chen JM, Fernandes R, Deering DW, Conley A (2005). “Methodology comparison for canopy structure parameters extraction from digital hemispherical photography in boreal forests.” Agricultural and Forest Meteorology, 129(3–4), 187–207. ISSN 0168-1923. doi:10.1016/j.agrformet.2004.09.006.

Macfarlane C (2011). “Classification method of mixed pixels does not affect canopy metrics from digital images of forest overstorey.” Agricultural and Forest Meteorology, 151(7), 833–840. doi:10.1016/j.agrformet.2011.01.019.

Ridler TW, Calvard S (1978). “Picture thresholding using an iterative selection method.” IEEE Transactions on Systems, Man, and Cybernetics, 8(8), 630–632. doi:10.1109/tsmc.1978.4310039.

Examples

## Not run: 
path <- system.file("external/DSCN4500.JPG", package = "rcaiman")
zenith_colrow <- c(1276, 980)
diameter <- 756*2
caim <- read_caim(path, zenith_colrow - diameter/2, diameter, diameter)
z <- zenith_image(ncol(caim), lens("Nikon_FCE9"))
r <- invert_gamma_correction(caim$Blue)
r <- correct_vignetting(r, z, c(0.0638, -0.101)) %>% normalize_minmax()
rings <- ring_segmentation(z, 15)
bin <- binarize_by_region(r, rings, "thr_twocorner")
plot(bin)

## End(Not run)

Description

Apply a threshold or a raster of thresholds to a grayscale image, producing a binary image.

Usage

binarize_with_thr(r, thr)

Arguments

Details

This function supports both global and pixel-wise thresholding. It is a wrapper around the > operator from the terra package. If a single numeric threshold is provided via thr, it is applied globally to all pixels in r. If instead a terra::SpatRaster object is provided, local thresholding is performed, where each pixel is compared to its corresponding threshold value.

This is useful after estimating thresholds using thr_twocorner(), thr_isodata(), or apply_by_direction(method = "thr_isodata"), among other posibilities.

Value

Logical terra::SpatRaster (TRUE for sky, FALSE for non-sky) with the same dimensions as r.

Note

For global thresholding, thr must be greater than or equal to the minimum value of r and lower than its maximum value.

Examples

r <- read_caim()
bin <- binarize_with_thr(r$Blue, thr_isodata(r$Blue[]))
plot(bin)

## Not run: 
# This function is also compatible with thresholds estimated using
# the 'autothresholdr' package:
require(autothresholdr)
r <- r$Blue
r <- normalize_minmax(r) %>% multiply_by(255) %>% round()
thr <- auto_thresh(r[], "IsoData")[1]
bin <- binarize_with_thr(r, thr)
plot(bin)

## End(Not run)

Description

Calculate the diameter in pixels of a 180 deg fisheye image.

Usage

calc_diameter(lens_coef, radius, angle)

Arguments

Details

This function is useful when the recording device has a field of view smaller than 180 deg. Given a lens projection function and data points consisting of radii (pixels) and their corresponding zenith angles ($\theta$), it returns the horizon radius (i.e., the radius for $\theta$ equal to 90 deg).

When working with non-circular hemispherical photography, this function helps determine the diameter that a circular image would have if the equipment recorded the whole hemisphere, required to build the correct zenith image to use as input for expand_noncircular().

The required data (radius–angle pairs) can be obtained following the instructions in the user manual of Hemisfer software. A slightly simpler alternative is:

1. Find a vertical wall and a leveled floor, both well-constructed.

2. Draw a triangle of $5 \times 4 \times 3$ meters on the floor, with the 4-meter side along the wall.

3. Place the camera over the vertex 3 meters away from the wall, at a chosen height (e.g., 1.3 m).

4. Make a mark on the wall at the chosen height over the wall-vertex nearest to the camera vertex. Make four more marks at 1 m intervals along a horizontal line. This creates marks for 0, 18, 34, 45, and 54 deg $\theta$.

5. Before taking the photograph, align the zenith coordinates with the 0 deg $\theta$ mark and ensure the optical axis is level.

The line selection tool of ImageJ can be used to measure the distance in pixels between points on the image. Draw a line and use the menu Analyze > Measure to obtain its length.

For obtaining the projection of a new lens, see calibrate_lens().

Value

Numeric vector of length one. Estimated diameter in pixels, rounded to the nearest even integer (see zenith_image() for details).

Examples

# Nikon D50 and Fisheye Nikkor 10.5mm lens
calc_diameter(lens("Nikkor_10.5mm"), 1202, 54)

Description

Convert zenith angles (degrees) to normalized radial distance using the lens projection model.

Usage

calc_relative_radius(angle, lens_coef)

Arguments

Details

This helper maps zenith angle(s) to a relative radius in [0, 1] given the lens projection coefficients.

Value

Numeric vector of the same length as angle, constrained to [0, 1].

Examples

calc_relative_radius(45, lens())

Description

Computes the angular distance, in radians, between directions defined by zenith and azimuth angles on the unit sphere.

Usage

calc_spherical_distance(z1, a1, z2, a2)

Arguments

Details

This function calculates the angle between two directions originating from the center of a unit sphere, using spherical trigonometry. The result is commonly referred to as spherical distance or angular distance. These terms are interchangeable when the sphere has radius one, as is standard in many applications, including celestial coordinate systems and, by extension, canopy hemispherical photography.

Spherical distance corresponds to the arc length of the shortest path between two points on the surface of a sphere. When the radius is one, this arc length equals the angle itself, expressed in radians.

Value

Numeric vector of the same length as z1 and a1, containing the spherical distance (in radians) from each (z1, a1) point to the reference direction (z2, a2).

Examples

set.seed(1)
z1 <- rnorm(10, 45, 20) * pi/180
a1 <- rnorm(10, 180, 90) * pi/180
calc_spherical_distance(z1, a1, 0, 0)

Description

Calculate zenith raster coordinates from points digitized with the open-source software package ‘ImageJ’.

Usage

calc_zenith_colrow(path_to_csv)

Arguments

Details

In this context, “zenith” denotes the location in the image that corresponds to the projection of the vertical direction when the optical axis is aligned vertically.

The technique described under the headline ‘Optical center characterization’ of the user manual of the software Can-Eye can be used to acquire the data for determining the zenith coordinates. This technique was used by Pekin and Macfarlane (2009), among others. Briefly, it consists in drilling a small hole in the cap of the fisheye lens (away from the center), and taking about ten photographs without removing the cap. The cap must be rotated about 30º before taking each photograph.

The point selection tool of ‘ImageJ’ software should be used to manually digitize the white dots and create a CSV file to feed this function. After digitizing the points on the image, use the dropdown menu Analyze>Measure to open the Results window. To obtain the CSV file, use File>Save As...

Another method (only valid when enough of the circle perimeter is depicted in the image) is taking a very bright picture (e.g., of a white-painted corner of a room) with the lens uncovered (do not use any mount). Then, digitize points over the circle perimeter. This was the method used for producing the example file (see Examples). It is worth noting that the perimeter of the circle depicted in a circular hemispherical photograph is not necessarily the horizon.

Value

Numeric vector of length two. Raster coordinates of the zenith. These coordinates follow image (raster) convention: the origin is in the upper-left, and the vertical axis increases downward, like a spreadsheet. This contrasts with Cartesian coordinates, where the vertical axis increases upward.

Note

This function assumes that all data points belong to the same circle, meaning that it does not support multiple holes when the Can-Eye procedure of drilling the lens cap is applied. The circle is fitted using the method presented by Kasa (1976).

References

Kasa I (1976). “A circle fitting procedure and its error analysis.” IEEE Transactions on Instrumentation and Measurement, IM–25(1), 8–14. ISSN 1557-9662. doi:10.1109/tim.1976.6312298.

Pekin B, Macfarlane C (2009). “Measurement of crown cover and leaf area index using digital cover photography and its application to remote sensing.” Remote Sensing, 1(4), 1298–1320. doi:10.3390/rs1041298.

Examples

## Not run: 
path <- system.file("external/points_over_perimeter.csv",
                    package = "rcaiman")
calc_zenith_colrow(path)

## End(Not run)

Description

Calibrate a fisheye lens to derive the mathematical relationship between image-space radial distances from the zenith and zenith angles in hemispherical space (assuming upward-looking hemispherical photography with the optical axis vertically aligned).

Usage

calibrate_lens(path_to_csv, degree = 3)

Arguments

Details

Fisheye lenses have a wide field of view and radial symmetry with respect to distortion. This property allows precise fitting of a polynomial model to relate pixel distances to zenith angles. The method implemented here, known as the "simple method", is described in detail by Díaz et al. (2024).

Value

List with named elements:

ds

Data frame used to fit the model.

model

lm object fitted to pixel distance vs. zenith angle.

horizon_radius

Radius at 90 deg.

lens_coef

Numeric vector of polynomial model coefficients for predicting relative radius.

zenith_colrow

Raster coordinates of the zenith push pin.

max_theta

Maximum zenith angle (deg).

max_theta_px

Distance in pixels between the zenith and the maximum zenith angle.

Step-by-step guide for producing a CSV file to feed this function

Materials

• this package and ImageJ

• camera and lens

• tripod

• standard yoga mat

• table at least as wide as the yoga mat width

• twenty two push pins of different colors

• one print of this sheet (A1 size, almost like a research poster).

• scissors

• some patience

Instructions

Cut the sheet by the dashed line. Place the yoga mat extended on top of the table. Place the sheet on top of the yoga mat. Align the dashed line with the yoga mat border closest to you. Place push pins on each cross. If you are gentle, the yoga mat will allow you to do that without damaging the table. Of course, other materials could be used to obtain the same result, such as cardboard, foam, nails, etc.

Place the camera on the tripod. Align its optical axis with the table while looking for getting an image showing the overlapping of the three pairs of push pins, as instructed in the print. In order to take care of the line of pins at 90º relative to the optical axis, it may be of help to use the naked eye to align the entrance pupil of the lens with the pins. The alignment of the push pins only guarantees the position of the lens entrance pupil, the leveling should be cheeked with an instrument, and the alignment between the optical axis and the radius of the zenith push pin should be taken into account. In practice, the latter is achieved by aligning the camera body with the orthogonal frame made by the quarter circle.

Take a photo and transfer it to the computer, open it with ImageJ, and use the point selection tool to digitize the push pins, starting from the zenith push pin and not skipping any shown push pin. End with an additional point where the image meets the surrounding black (or the last pixel in case there is not blackness because it is not a circular hemispherical image. There is no need to follow the line formed by the push pins). Then, use the dropdown menu Analyze>Measure to open the window Results. To obtain the CSV, use File>Save As...

Note

To calibrate different directions, think of the fisheye image as an analog clock. To calibrate 3 o'clock, attach the camera to the tripod in landscape mode while leaving the quarter-circle at the lens's right side. To calibrate 9 o'clock, rotate the camera to put the quarter-circle at the lens's left side. To calibrate 12 and 6 o'clock, do the same but with the camera in portrait mode.

References

Díaz GM, Lang M, Kaha M (2024). “Simple calibration of fisheye lenses for hemipherical photography of the forest canopy.” Agricultural and Forest Meteorology, 352, 110020. ISSN 0168-1923. doi:10.1016/j.agrformet.2024.110020.

See Also

test_lens_coef(), crosscalibrate_lens(), extract_radiometry()

Examples

path <- system.file("external/Results_calibration.csv", package = "rcaiman")
calibration <- calibrate_lens(path)
coefficients(calibration$model)
calibration$lens_coef %>% signif(3)
calibration$horizon_radius

## Not run: 
test_lens_coef(calibration$lens_coef) #MacOS and Windows tend to differ here
test_lens_coef(c(0.628, 0.0399, -0.0217))

## End(Not run)

.fp <- function(theta, lens_coef) {
  x <- lens_coef[1:5]
  x[is.na(x)] <- 0
  for (i in 1:5) assign(letters[i], x[i])
  a * theta + b * theta^2 + c * theta^3 + d * theta^4 + e * theta^5
}

plot(calibration$ds)
theta <- seq(0, pi/2, pi/180)
lines(theta, .fp(theta, coefficients(calibration$model)))

Description

Segment a raster into square regions of equal size arranged in a chessboard-like pattern.

Usage

chessboard(r, size)

Arguments

Details

This function divides the extent of r into non-overlapping square segments of the given size, producing a segmentation map where each segment has a unique integer label. It can be an alternative to skygrid_segmentation() in special cases.

Value

terra::SpatRaster with one layer and integer values, where each unique value corresponds to a square-segment ID.

Examples

caim <- read_caim()
seg <- chessboard(caim, 20)
plot(caim$Blue)
plot(extract_feature(caim$Blue, seg))

Description

Generate an image of relative radiance or luminance based on the CIE General Sky model.

Usage

cie_image(z, a, sun_angles, sky_coef)

Arguments

Value

terra::SpatRaster with one layer whose pixel values represent relative luminance or radiance across the sky hemisphere, depending on whether the data used to obtain sky_coef was luminance or radiance.

Note

Coefficient sets and formulation are available in cie_table.

Examples

z <- zenith_image(50, lens())
a <- azimuth_image(z)
sky_coef <- cie_table[4,1:5] %>% as.numeric()
sun_angles <- c(z = 45, a = 0)
plot(cie_image(z, a, sun_angles, sky_coef))

Description

The Commission Internationale de l’Éclairage (CIE; International Commission on Illumination) standard (CIE 2004) defines 15 pre-calibrated sky luminance distributions, each described by a pair of analytical functions, the gradation function $\Phi(\theta) = 1 + a \cdot \exp\left(\frac{b}{\cos\theta}\right)$, and the indicatrix function $f(\chi) = 1 + c \cdot \left[ \exp(d \cdot \chi) - \exp\left(d \cdot \frac{\pi}{2}\right) \right] + e \cdot \cos^2\chi$. Combined, they can predict the relative radiance $\rho_R$ in any sky direction $(\theta, \phi)$ as: $\hat{\rho_R}^{\circ}(\theta, \phi) = \frac{\Phi(\theta) \cdot f(\chi(\theta, \phi; \theta_\odot, \phi_\odot))} {\Phi(0) \cdot f(\chi(\theta, \phi; 0, 0))}$, where $\theta$ is the zenith angle, $\phi$ is the azimuth angle, and $\theta_\odot, \phi_\odot$ are the zenith and azimuth of the sun disk.

Usage

cie_table

Format

data.frame with 15 rows and 8 columns:

a

gradation function parameter.

b

gradation function parameter.

c

indicatrix function parameter.

d

indicatrix function parameter.

e

indicatrix function parameter.

indicatrix_group

factor with six categories and numerical tags.

general_sky_type

factor with three categories: "Overcast", "Clear", and "Partly cloudy".

description

user-friendly description of the sky type.

Source

Li et al. (2016)

References

CIE (2004). “ISO 15469:2004(E) / CIE S 011/E:2003 - Spatial distribution of daylight — CIE standard general sky.” https://www.iso.org/standard/38608.html. International Standard.

Li DHW, Lou S, Lam JC, Wu RHT (2016). “Determining solar irradiance on inclined planes from classified CIE (International Commission on Illumination) standard skies.” Energy, 101, 462–470. doi:10.1016/j.energy.2016.02.054.

Description

Compute three color-opponent gradients to enhance the visual separation between sky and canopy in hemispherical photographs, particularly under diffuse light or complex cloud patterns.

Usage

complementary_gradients(caim)

Arguments

Details

The method exploits chromatic differences between the red, green, and blue bands, following a simplified opponent-color logic. Each gradient is normalized by total brightness and modulated by a logistic contrast function to reduce the influence of underexposed regions:

• "green_magenta" = $(R - G + B) / (R + G + B)$ · logistic(brightness)

• "yellow_blue" = $(-R - G + B) / (R + G + B)$ · logistic(brightness)

• "red_cyan" = $(-R + G + B) / (R + G + B)$ · logistic(brightness)

The logistic(brightness) term is computed as:

$\text{logistic}(x) = \frac{1}{1 + \exp\left(-\frac{x - q_{0.1}}{\mathrm{IQR}}\right)}$

where $q_{0.1}$ is the 10th percentile of brightness values ($x = R + G + B$), and $IQR$ is their interquartile range.

This weighting suppresses gradients in poorly exposed regions to reduce spurious values caused by low signal-to-noise ratios.

Value

Numeric terra::SpatRaster with three layers and the same geometry as caim. The layers ("green_magenta", "yellow_blue", "red_cyan") are chromatic gradients modulated by brightness.

Note

This function is part of a paper under preparation.

References

There are no references for Rd macro ⁠\insertAllCites⁠ on this help page.

Examples

## Not run: 
caim <- read_caim()
com <- complementary_gradients(caim)
plot(com)

## End(Not run)

Description

Calculate canopy openness from a binarized hemispherical image with angular coordinates.

Usage

compute_canopy_openness(bin, z, a, m = NULL, angle_width = 10)

Arguments

Details

Canopy openness is computed following the equation from Gonsamo et al. (2011):

$CO = \sum_{i = 1}^{N} GF(\phi_i, \theta_i) \cdot \frac{\cos(\theta_{1,i}) - \cos(\theta_{2,i})}{n_i}$

where $GF(\phi_i, \theta_i)$ is the gap fraction in cell $i$, $\theta_{1,i}$ and $\theta_{2,i}$ are the lower and upper zenith angles of the cell, $n_i$ is the number of cells in the corresponding zenith ring, and $N$ is the total number of cells.

When a mask is provided via the m argument, the equation is adjusted to compensate for the reduced area of the sky vault:

$CO = \frac{\sum_{i=1}^{N} GF(\phi_i, \theta_i) \cdot w_i} {\sum_{i=1}^{N} w_i} \quad \text{with} \quad w_i = \frac{\cos(\theta_{1,i}) - \cos(\theta_{2,i})}{n_i}$

The denominator ensures that the resulting openness value remains scale-independent. Without this normalization, masking would lead to underestimation, as the numerator alone assumes full hemispherical coverage.

Value

Numeric vector of length one, constrained to the range $[0, 1]$.

References

Gonsamo A, Walter JN, Pellikka P (2011). “CIMES: A package of programs for determining canopy geometry and solar radiation regimes through hemispherical photographs.” Computers and Electronics in Agriculture, 79(2), 207–215. doi:10.1016/j.compag.2011.10.001.

Examples

caim <- read_caim()
z <- zenith_image(ncol(caim), lens())
a <- azimuth_image(z)
m <- select_sky_region(z, 0, 70)
bin <- binarize_with_thr(caim$Blue, thr_isodata(caim$Blue[m]))
plot(bin)
compute_canopy_openness(bin, z, a, m, 10)

Description

Create an RGB image that resembles a photo taken with a conventional lens, using a small patch from the example hemispherical image.

Usage

conventional_lens_image()

Details

This is a fixed crop and reorientation of read_caim(). It does not perform any re-projection. Intended for documentating functions.

The following code was used to define the region:

caim <- read_caim()
r <- caim$Blue
z <- zenith_image(ncol(caim), lens())
a <- azimuth_image(z)

m <- rast(z)
m[] <- calc_spherical_distance(
  z[] * pi / 180,
  a[] * pi / 180,
  1,# hinge-angle
  90 * pi / 180
)
m <- !binarize_with_thr(m, 30 * pi / 180)
m[is.na(z)] <- 0
m

x11()
plot(m * caim$Blue)
za <- click(c(z, a))
za
row_col <- row_col_from_zenith_azimuth(z, a, za[,1], za[,2])
plot(caim$Blue)
points(row_col$col, nrow(caim) - row_col$row, col = 2, pch = 10)
mn_y <- min(nrow(caim) -row_col$row)
mx_y <- max(nrow(caim) -row_col$row)
mn_x <- min(row_col$col)
mx_x <- max(row_col$col)
r <- terra::crop(caim$Blue, terra::ext(mn_x, mx_x, mn_y, mx_y))
plot(r)

Value

Three-layer terra::SpatRaster with bands in RGB order.

See Also

read_caim()

Examples

conventional_lens_image()

Description

Apply a vignetting correction to an image using a polynomial model.

Usage

correct_vignetting(
  r,
  z,
  lens_coef_v,
  method = "simple",
  model_type = "polynomial"
)

Arguments

Details

Vignetting is the gradual reduction of image brightness toward the periphery. This function applies a parametric correction tailored to a specific hemispherical camera system that varies with zenith angle at the pixel level.

The selected calibration scheme defines the modeling conventions used by rcaiman to represent and correct vignetting effects. Different schemes correspond to different assumptions and degrees of freedom in the underlying polynomial model.

simple

A constrained polynomial scheme in which the vignetting function is normalized such that $f_v(0) = 1$. The effective model is $f_v(\theta) = 1 + a\theta + b\theta^2 + \dots + m\theta^n$, where $\theta$ is the zenith angle (in radians) and $a,b,\dots,m$ are the polynomial coefficients. Polynomial degrees up to 6 are supported. This scheme follows the approach described in Díaz et al. (2024). See extract_radiometry() for guidance on estimating these coefficients.

free_form

A flexible polynomial scheme in which all model parameters, including the constant term, are estimated from data. The effective model is $f_v(\theta) = a + b\theta + c\theta^2 + \dots + m\theta^n$, where $\theta$ is the zenith angle (in radians) and $a,b,\dots,m$ are free parameters. Up to 7 parameters are supported. This scheme assumes that the input data reliably capture the true shape of the vignetting function, and therefore applies no structural constraints on the model parameters. It is typically used with data acquired using a photometric sphere, following the methodology described in Lang et al. (2010).

Value

terra::SpatRaster with the same content as r but with pixel values adjusted to correct for vignetting, preserving all other properties (layers, names, extent, and CRS).

References

Díaz GM, Lang M, Kaha M (2024). “Simple calibration of fisheye lenses for hemipherical photography of the forest canopy.” Agricultural and Forest Meteorology, 352, 110020. ISSN 0168-1923. doi:10.1016/j.agrformet.2024.110020.

Lang M, Kuusk A, Mõttus M, Rautiainen M, Nilson T (2010). “Canopy gap fraction estimation from digital hemispherical images using sky radiance models and a linear conversion method.” Agricultural and Forest Meteorology, 150(1), 20–29. doi:10.1016/j.agrformet.2009.08.001.

Examples

## Not run: 
path <- system.file("external/APC_0836.jpg", package = "rcaiman")
caim <- read_caim(path)
z <- zenith_image(2132, lens("Olloclip"))
a <- azimuth_image(z)
zenith_colrow <- c(1063, 771)

caim <- expand_noncircular(caim, z, zenith_colrow)
m <- !is.na(caim$Red) & !is.na(z)
caim[!m] <- 0

bin <- binarize_with_thr(caim$Blue, thr_isodata(caim$Blue[m]))
display_caim(caim$Blue, bin)

caim <- invert_gamma_correction(caim, 2.2)
caim <- correct_vignetting(caim, z, c(-0.0546, -0.561, 0.22)) %>%
        normalize_minmax()

## End(Not run)

Description

Apply vignetting correction using the aperture-dependent empirical model described by Lang et al. (2010).

Usage

correct_vignetting_legacy(r, z, aperture)

Arguments

Details

This function implements the aperture-dependent vignetting model reported in Eq. (5) of Lang et al. (2010). The relative pixel value with respect to the image center is given by:

$B = 0.622 - \frac{1.504 \theta_v}{A^{4.553}} + 0.0539 \ln\left((22.87 - A)(91.0 - \theta_v)\right)$

where $\theta_v$ is the zenith angle in degrees (0–90 deg) and $A$ is the aperture value used during image acquisition.

This function is provided to support workflows based on that historical calibration. It does not use the polynomial schemes implemented in correct_vignetting().

Value

Numeric terra::SpatRaster with the same geometry as r, with pixel values divided by $B$ to correct for vignetting.

References

Lang M, Kuusk A, Mõttus M, Rautiainen M, Nilson T (2010). “Canopy gap fraction estimation from digital hemispherical images using sky radiance models and a linear conversion method.” Agricultural and Forest Meteorology, 150(1), 20–29. doi:10.1016/j.agrformet.2009.08.001.

Description

Extracts a rectangular region of interest (ROI) from a canopy image. This function complements read_caim() and read_caim_raw().

Usage

crop_caim(r, upper_left = NULL, width = NULL, height = NULL)

Arguments

Value

terra::SpatRaster object containing the same layers and values as r but restricted to the selected ROI, preserving all other properties.

Selecting a Region of Interest

To load a specific subregion from the image, use the arguments upper_left, width, and height. These are expressed in raster coordinates, similar to a spreadsheet layout: columns first, then rows. In other words, specify coordinates as c(column, row), not c(row, column), which is typical in data.frame objects.

While any image editor can be used to obtain these values, this function was tested with ImageJ, particularly the Fiji distribution. A recommended workflow:

1. Open the image in Fiji.

2. Draw a rectangular selection.

3. Go to Edit > Selection > Specify... to read upper_left, width, and height.

Note

rcaiman uses terra without geographic semantics: rasters are kept with unit resolution (cell size = 1) and a standardized extent ext(0, ncol, 0, nrow) with CRS EPSG:7589 (set_rcaiman_geometry()).

Examples

caim <- read_caim()
ncell(caim)
caim <- crop_caim(caim, c(231,334), 15, 10)
ncell(caim)

path <- system.file("external/DSCN4500.JPG", package = "rcaiman")
caim <- read_caim(path)
zenith_colrow <- c(1286, 986)
horizon_radius <- 742
upper_left <- zenith_colrow - horizon_radius
height <- width <- horizon_radius * 2
plot(crop_caim(caim, upper_left, width, height))
display_caim(crop_caim(caim, upper_left, width, height))

Description

Given two photographs taken from the same point (matching entrance pupils and aligned optical axes), with calibrated and uncalibrated cameras, derives a polynomial projection for the uncalibrated device. Intended for cases where a camera calibrated with a method of higher accuracy than calibrate_lens() is available, or when there is a main camera to which all other devices should be adjusted.

Points must be digitized in tandem with ImageJ and saved as CSV files. See calibrate_lens() for background and general concepts.

Usage

crosscalibrate_lens(
  path_to_csv_uncal,
  path_to_csv_cal,
  zenith_colrow_uncal,
  zenith_colrow_cal,
  diameter_cal,
  lens_coef,
  degree = 3
)

Arguments

Details

Estimate a lens projection for an uncalibrated camera by referencing a calibrated camera photographed from the exact same location.

Value

List with components:

ds

data.frame with zenith angle (theta, radians) and pixel radius (px) from the uncalibrated camera.

model

lm object: polynomial fit of px ~ theta.

horizon_radius

numeric vector of length one. Pixel radius at 90 deg.

lens_coef

numeric vector. Distortion coefficients normalized by horizon_radius.

See Also

calibrate_lens(), calc_zenith_colrow()

Description

Converts fuzzy membership values into a binary classification using a regional approach that preserves aggregation consistency between the fuzzy and binary representations.

Usage

defuzzify(mem, seg)

Arguments

Details

The conversion is applied within segments defined by seg, ensuring that, in each segment, the aggregated Boolean result matches the aggregated fuzzy value. This approach is well suited for converting subpixel estimates, such as gap fraction, into binary outputs.

Value

Logical terra::SpatRaster of the same dimensions as mem, where each pixel value represents the binary version of mem after applying the regional defuzzification procedure.

Note

This method is also available in the HSP software package. See hsp_compat().

Examples

## Not run: 
caim <- read_caim()
r <- caim$Blue
z <- zenith_image(ncol(caim), lens())
a <- azimuth_image(z)

path <- system.file("external/example.txt", package = "rcaiman")
sky_cie <- read_sky_cie(gsub(".txt", "", path), z, a)

sky_above <- ootb_sky_above(sky_cie$model$rr$sky_points, z, a, sky_cie)

ratio <- r / sky_above$dn_raster
ratio <- normalize_minmax(ratio, 0, 1, TRUE)
plot(ratio)
g <- skygrid_segmentation(z, a, 10)
bin2 <- defuzzify(ratio, g)
plot(bin2) # unsatisfactory results due to light conditions

## End(Not run)

Description

Wrapper for EBImage::display() that streamlines the visualization of canopy images, optionally overlaying binary masks and segmentation borders. It is intended for quick inspection of processed or intermediate results in a graphical viewer.

Usage

display_caim(
  caim = NULL,
  bin = NULL,
  seg = NULL,
  sampling_points = NULL,
  sun_row_col = NULL,
  sun_disk_size = 9
)

Arguments

Value

Invisible NULL. Called for side effects (image viewer popup).

Examples

## Not run: 
caim <- read_caim()
r <- caim$Blue
z <- zenith_image(ncol(caim), lens())
a <- azimuth_image(z)

# See fit_cie_model() for details on below file
path <- system.file("external/sky_points.csv",
                    package = "rcaiman")
sky_points <- read.csv(path)
sky_points <- sky_points[c("Y", "X")]
colnames(sky_points) <- c("row", "col")
head(sky_points)
plot(caim$Blue)
points(sky_points$col, nrow(caim) - sky_points$row, col = 2, pch = 10)

x11()
# plot(caim$Blue)
# sun_angles <- click(c(z, a), 1) %>% as.numeric()
sun_angles <- c(z = 49.5, a = 27.4) #taken with above lines then hardcoded

sun_row_col <- row_col_from_zenith_azimuth(z, a,
                                           sun_angles["z"],
                                           sun_angles["a"])

bin <- binarize_with_thr(r$Blue, thr_isodata(r$Blue[]))
seg <- equalarea_segmentation(z, a, 200)

# color image
display_caim(caim)
# greyscale
display_caim(caim$Blue)
# to check binarization quality (press `h` with the pointer on the image to
# learn useful hotkeys)
display_caim(caim$Blue, bin)
# to see the segments
display_caim(seg = seg)
# to see the marks
display_caim(caim, sampling_points = sky_points, sun_row_col = sun_row_col)

## End(Not run)

Description

Segment a hemispherical view into n_cells equal-area cells in spherical space.

Usage

equalarea_segmentation(z, a, n_cells)

Arguments

Details

Segmentation begins by creating n_cells initial rings through uniform sampling in cos(z) space. The core logic of the algorithm is grouping k initial rings into a larger ring and then subdividing that ring azimuthally into k cells. This procedure ensures equal area while allowing the final segments to adopt compact shapes. In this way, the larger rings become the final rings, of which there will be n_rings = round(sqrt(n_cells)).

The final ring containing the zenith has k = 1, the ring touching the horizon has k = n_rings, and the intermediate k values must be determined. To impose a compactness structure across zenith angles, the method uses the equisolid-angle projection:

$r = 2 \, f \, \sin\!\left(\frac{\theta}{2}\right)$

where $r$ is the radial distance, $f$ is the focal length, and $\theta$ is the zenith angle. Based on this transformation, the expression round(2 * sin(theta / 2) * n_rings * x) is used to compute k, with x being a scalar (estimated with optim()) that adjusts the sum of k values to exactly match n_cells.

Because the true central zenith angles of the final rings are not known in advance (because they depend on k) the procedure is iterated. It starts by assuming rings of equal angular width to obtain provisional central zenith angles, computes k, updates the ring boundaries from the resulting grouping, recomputes the central zenith angles, and repeats until convergence.

Value

Single-layer terra::SpatRaster with integer values from 1 to n_cells.

See Also

skygrid_segmentation()

Examples

## Not run: 
z <- zenith_image(500, lens())
a <- azimuth_image(500, lens())
seg <- equalarea_segmentation(z, a, n_cells = 512)
plot(seg)

## End(Not run)

Description

Estimates the sun’s zenith and azimuth angles (deg) from a canopy hemispherical photograph, using either direct detection of the solar disk or indirect cues from the circumsolar region.

Usage

estimate_sun_angles(
  r,
  z,
  a,
  bin,
  seg,
  angular_radius_sun = 30,
  mode = "obscured"
)

Arguments

Details

This function can operate under two alternative assumptions for estimating the sun position:

Veiled sun

The solar disk is visible or partially obscured; the position is inferred from localized brightness peaks. It is also useful when a sunscreen is used, as in Kuusk and Paas (2007).

Obscured sun

The solar disk is not visible; the position is inferred from radiometric and spatial cues aggregated over the circumsolar region.

When mode = "veiled", seg and angular_radius_sun are ignored. Estimates refer to positions above the horizon; therefore, estimated angles may require further manipulation if the photograph was acquired under crepuscular light.

Value

Named numeric vector of length two, with names z and a, representing the sun’s zenith and azimuth angles (in degrees).

Note

A scientific article presenting and validating this method is currently under preparation.

References

Kuusk A, Paas M (2007). “Radiometric correction of hemispherical images.” ISPRS Journal of Photogrammetry and Remote Sensing, 61(6), 405–413. ISSN 0924-2716. doi:10.1016/j.isprsjprs.2006.10.005.

Examples

## Not run: 
caim <- read_caim()
z <- zenith_image(ncol(caim), lens())
a <- azimuth_image(z)
m <- !is.na(z)
r <- caim$Blue

bin <- binarize_by_region(r, ring_segmentation(z, 15), "thr_isodata") &
  select_sky_region(z, 0, 88)

seg <- skygrid_segmentation(z, a, 10)
sun_angles <- estimate_sun_angles(r, z, a, bin, seg,
                                  angular_radius_sun = 30)
row_col <- row_col_from_zenith_azimuth(z, a,
                                       sun_angles["z"],
                                       sun_angles["a"])
plot(caim$Blue)
points(row_col[1,2], nrow(caim) - row_col[1,1], col = "yellow",
       pch = 8, cex = 3)

## End(Not run)

Description

Add NA margins to a hemispherical photograph to align radiance at the zenith with the image center. In this context, “zenith” denotes the location in the image that corresponds to the projection of the vertical direction when the optical axis is aligned vertically. Intended for non-circular images.

Usage

expand_noncircular(caim, z, zenith_colrow)

Arguments

Value

terra::SpatRaster with the same layers and pixel values as caim, but with NA margins added to center the zenith.

Note

rcaiman uses terra without geographic semantics: rasters are kept with unit resolution (cell size = 1) and a standardized extent ext(0, ncol, 0, nrow) with CRS EPSG:7589.

Examples

## Not run: 
# Non-circular fisheye images from a smartphone with an auxiliary Lens
# (also applicable to non-circular fisheye images from DSLR cameras)
path <- system.file("external/APC_0836.jpg", package = "rcaiman")
caim <- read_caim(path)
z <- zenith_image(2132,  lens("Olloclip"))
a <- azimuth_image(z)
zenith_colrow <- c(1063, 771)
caim <- expand_noncircular(caim, z, zenith_colrow)
plot(caim$Blue, col = seq(0, 1, 1/255) %>% grey())
m <- !is.na(caim$Red) & !is.na(z)
plot(m, add = TRUE, alpha = 0.3, legend = FALSE)

## End(Not run)

Description

Compute a robust coefficient of variation (CV) at each sky point using a $3 \times 3$ neighbourhood in a raster image.

Usage

extract_cv(r, sky_points)

Arguments

Details

For each sky point, the function extracts the radiance values in a $3 \times 3$ window centered on the corresponding pixel of r, and computes a robust CV defined as:

$\mathrm{CV}_\mathrm{robust} = \mathrm{MAD} / \mathrm{median}$

where MAD is the median absolute deviation. This measure is suited for detecting local fluctuations caused by bright canopy elements or mixed pixels, while being less sensitive to outliers than a mean–sd based CV.

The extraction procedure mirrors the spatial logic used in extract_dn(), but instead of averaging the neighbourhood, it applies a robust dispersion measure (MAD/median) to quantify local radiance variability around each point.

Value

Numeric vector of length equal to nrow(sky_points), containing one robust CV value per sky point.

See Also

extract_dn() for extracting digital numbers and tune_sky_sampling() for an application of this metric.

Examples

## Not run: 
caim <- read_caim()
r <- caim$Blue
z <- zenith_image(ncol(caim), lens())
a <- azimuth_image(z)

# See fit_cie_model() for details on below file
path <- system.file("external/sky_points.csv",
                    package = "rcaiman")
sky_points <- read.csv(path)
sky_points <- sky_points[c("Y", "X")]
colnames(sky_points) <- c("row", "col")
head(sky_points)
plot(caim$Blue)
points(sky_points$col, nrow(caim) - sky_points$row, col = 2, pch = 10)

extract_cv(caim$Blue, sky_points)


## End(Not run)

Description

Obtain digital numbers from a raster at positions defined by points, with optional local averaging.

Usage

extract_dn(r, sampling_points, use_window = TRUE)

Arguments

Details

Wraps terra::extract() to support a $3 \times 3$ window centered on each target pixel (local mean). When it is disabled, only the central pixel value is retrieved.

Value

data.frame containing the original sampling_points plus one column per layer in r (named after the layers).

Note

For instructions on manually digitizing points, see the “Digitizing sky points with ImageJ” and “Digitizing sky points with QGIS” sections in fit_cie_model().

See Also

sample_sky_points() for automatically sampling sky points.

Examples

## Not run: 
caim <- read_caim()
r <- caim$Blue
z <- zenith_image(ncol(caim), lens())
a <- azimuth_image(z)

# See fit_cie_model() for details on below file
path <- system.file("external/sky_points.csv",
                    package = "rcaiman")
sky_points <- read.csv(path)
sky_points <- sky_points[c("Y", "X")]
colnames(sky_points) <- c("row", "col")
head(sky_points)
plot(caim$Blue)
points(sky_points$col, nrow(caim) - sky_points$row, col = 2, pch = 10)

sky_points <- extract_dn(caim, sky_points)
head(sky_points)

# To aggregate DNs across points (excluding 'row' and 'col'):
apply(sky_points[, -(1:2)], 2, mean, na.rm = TRUE)

## End(Not run)

Description

Extract a numeric or logical summary from segmented raster regions using a user-defined reducer, returning one value per segment as a raster map or a named vector.

Usage

extract_feature(r, seg, fun = mean, return = "raster", ignore_label_0 = TRUE)

Arguments

Details

Segments labeled 0 can be ignored via ignore_label_0 = TRUE. The function in fun must return a single numeric or logical value for any input vector (e.g., mean, median, or a custom reducer).

Value

If return = "raster", a terra::SpatRaster where each pixel holds its segment’s feature value. If return = "vector", a named numeric (or logical) vector with one value per segment.

Examples

r <- read_caim()
z <- zenith_image(ncol(r),lens())
a <- azimuth_image(z)
g <- skygrid_segmentation(z, a, 10)
print(extract_feature(r$Blue, g, return = "vector"))
# plot(extract_feature(r$Blue, g, return = "raster"))

Description

Buid a datasets for vignetting modeling by sistematically extracting radiometry from images taken with the aid of a portable light source and the calibration board detailed in calibrate_lens().

Usage

extract_radiometry(l, size_px = NULL)

Arguments

Value

data.frame with two columns:

theta

zenith angle in radians.

radiometry

normalized digital number.

Guidance

Lenses have the inconvenient property of increasingly attenuating light along the direction orthogonal to the optical axis. This phenomenon is known as the vignetting effect and varies with lens model and aperture setting. The method outlined here, known as the simple method, is explained in details in Díaz et al. (2024). Next explanation might serve mostly as a quick recap of it.

The development of the simple method was done with a Kindle Paperwhite eBooks reader of 6" with built-in light. However, an iPhone 6 plus was also tested in the early stages of development and no substantial differences were observed. A metal bookends desk book holder was used to fasten the eBook reader upright and a semi-transparent paper to favor a Lambertian light distribution. In addition, the latter was used to draw on in order to guide pixel sampling. The book holder also facilitated the alignment of the screen with the dotted lines of the printed quarter-circle.

As a general guideline, a wide variety of mobile devices could be used as light sources, but if scattered data points are obtained with it, then other light sources should be tested in order to double check that the light quality is not the reason for points scattering.

With the room only illuminated by the portable light source, nine photographs should be taken with the light source located in the equivalent to 0, 10, 20, 30, 40, 50, 60, 70, and 80 degrees of zenith angle, respectively. Camera configuration should be in manual mode and set with the aperture (f/number) for which a vignetting function is required. The shutter speed should be regulated to obtain light-source pixels with middle grey values. The nine photographs should be taken without changing the camera configuration and the light conditions.

This code exemplifies how to use the function to obtain data and base R functions to obtain the vignetting function ($f_v$).

zenith_colrow <- c(1500, 997)
diameter <- 947*2
z <- zenith_image(diameter, c(0.689, 0.0131, -0.0295))
a <- azimuth_image(z)

.read_raw <- function(path_to_raw_file) {
  r <- read_caim_raw(path_to_raw_file, only_blue = TRUE)
  r <- crop_caim(r, zenith_colrow - diameter/2, diameter, diameter)
  r <- fisheye_to_equidistant(r, z, a, m, radius = diameter/2,
                              k = 1, p = 1, rmax = 100)
}

l <- Map(.read_raw, dir("raw/up/", full.names = TRUE))
up_data <- extract_radiometry(l)
l <- Map(.read_raw, dir("raw/down/", full.names = TRUE))
down_data <- extract_radiometry(l)
l <- Map(.read_raw, dir("raw/right/", full.names = TRUE))
right_data <- extract_radiometry(l)
l <- Map(.read_raw, dir("raw/left/", full.names = TRUE))
left_data <- extract_radiometry(l)

ds <- rbind(up_data, down_data, right_data, left_data)

plot(ds, xlim = c(0, pi/2), ylim= c(0.5,1.05),
      col = c(rep(1,9),rep(2,9),rep(3,9),rep(4,9)))
legend("bottomleft", legend = c("up", "down", "right", "left"),
       col = 1:4, pch = 1)

fit <- lm((1 - ds$radiometry) ~ poly(ds$theta, 3, raw = TRUE) - 1)
summary(fit)
coef <- -fit$coefficients #did you notice the minus sign?
.fv <- function(x) 1 + coef[1] * x + coef[2] * x^2 + coef[3] * x^3
curve(.fv, add = TRUE, col = 2)
coef

Once one of the aperture settings is calibrated, it can be used to calibrate all the rest. To do so, the equipment should be used to take photographs in all desired exposition and without moving the camera, including the aperture already calibrated and preferably under an open sky in stable diffuse light conditions. Below code can be used as a template.

zenith_colrow <- c(1500, 997)
diameter <- 947*2
z <- zenith_image(diameter, c(0.689, 0.0131, -0.0295))
a <- azimuth_image(z)

files <- dir("raw/", full.names = TRUE)
l <- list()
for (i in seq_along(files)) {
  if (i == 1) {
    # because the first aperture was the one already calibrated
    .read_raw <- function(path_to_raw_file) {
        r <- read_caim_raw(path_to_raw_file, only_blue = TRUE)
        r <- crop_caim(r, zenith_colrow - diameter/2, diameter, diameter)
        r <- correct_vignetting(r, z, c(0.0302, -0.320, 0.0908))
        r <- fisheye_to_equidistant(r, z, a, m, radius = diameter/2,
                                    k = 1, p = 1, rmax = 100)
      }
  } else {
      .read_raw <- function(path_to_raw_file) {
          r <- read_caim_raw(path_to_raw_file, only_blue = TRUE)
          r <- crop_caim(r, zenith_colrow - diameter/2, diameter, diameter)
          r <- fisheye_to_equidistant(r, z, a, m, radius = diameter/2,
                                      k = 1, p = 1, rmax = 100)
        }
  }
  l[[i]] <- .read_raw(files[i])
}

ref <- l[[1]]
rings <- ring_segmentation(zenith_image(ncol(ref), lens()), 3)
theta <- seq(1.5, 90 - 1.5, 3) * pi/180

.fun <- function(r) {
  r <- extract_feature(r, rings, return = "vector")
  r/r[1]
}

l <- Map(.fun, l)

.fun <- function(x) {
  x / l[[1]][] # because the first is the one already calibrated
}
radiometry <- Map(.fun, l)

l <- list()
for (i in 2:length(radiometry)) {
  l[[i-1]] <- data.frame(theta = theta, radiometry = radiometry[[i]][])
}
ds <- l[[1]]
head(ds)

The result is one dataset (ds) for each file. This is all what it is needed before using base R functions to fit a vignetting function

Note

This function does not fit the vignetting function itself. The output is a dataset to be used in subsequent modeling steps. See above sections for guidance.

References

Díaz GM, Lang M, Kaha M (2024). “Simple calibration of fisheye lenses for hemipherical photography of the forest canopy.” Agricultural and Forest Meteorology, 352, 110020. ISSN 0168-1923. doi:10.1016/j.agrformet.2024.110020.

Description

Compute relative radiance at selected sky points by dividing their digital numbers (DN) by an estimated zenith DN.

Usage

extract_rr(r, z, a, sky_points, no_of_points = 3, use_window = TRUE)

Arguments

Value

List with named elements:

zenith_dn

numeric. Estimated DN at the zenith.

sky_points

data.frame with columns row, col, a, z, dn, and rr (pixel location, angular coordinates, extracted DN, and relative radiance). If no_of_points is NULL, zenith_dn = 1 and dn = rr.

Examples

## Not run: 
caim <- read_caim()
z <- zenith_image(ncol(caim), lens())
a <- azimuth_image(z)

# See fit_cie_model() for details on the CSV file
path <- system.file("external/sky_points.csv",
                    package = "rcaiman")
sky_points <- read.csv(path)
sky_points <- sky_points[c("Y", "X")]
colnames(sky_points) <- c("row", "col")
head(sky_points)

plot(caim$Blue)
points(sky_points$col, nrow(caim) - sky_points$row, col = 2, pch = 10)
rr <- extract_rr(caim$Blue, z, a, sky_points, 1)
points(rr$sky_points$col, nrow(caim) - rr$sky_points$row, col = 3, pch = 0)

## End(Not run)

Description

Return image row and column indices for points distributed on the upper hemisphere using a spherical Fibonacci lattice with approximately constant angular spacing.

Usage

fibonacci_points(z, a, spacing)

Arguments

Value

data.frame with integer columns row and col.

See Also

skygrid_centers()

Examples

z <- zenith_image(100, lens())
a <- azimuth_image(z)
sampling_points <- fibonacci_points(z, a, 30)
## Not run: 
display_caim(is.na(z), sky_points = sky_points)

## End(Not run)

Description

Reproject a hemispherical image from fisheye to equidistant projection (also known as polar projection) to standardize its geometry for subsequent analysis and comparison between images.

Usage

fisheye_to_equidistant(r, z, a, m, radius = NULL, k = 1, p = 1, rmax = 10)

Arguments

Details

Pixel values and coordinates are treated as 3D points and reprojected using Cartesian interpolation. Internally, this function uses lidR::knnidw() as interpolation engine, so arguments k, p, and rmax are passed to it without modification.

Value

terra::SpatRaster with the same number of layers as r, reprojected to equidistant projection with circular shape and radius given by 'radius

Examples

## Not run: 
path <- system.file("external/APC_0836.jpg", package = "rcaiman")
caim <- read_caim(path)
calc_diameter(c(0.801, 0.178, -0.179), 1052, 86.2)
z <- zenith_image(2216,  c(0.801, 0.178, -0.179))
a <- azimuth_image(z)
zenith_colrow <- c(1063, 771)

caim <- expand_noncircular(caim, z, zenith_colrow)
caim2 <- fisheye_to_equidistant(caim, z, a, !is.na(z), radius = 600)
m <- !is.na(caim$Red) & select_sky_region(z, 0, 86.2)
caim[!m] <- 0
m2 <- fisheye_to_equidistant(m, z, a, !is.na(z), radius = 600)
m2 <- binarize_with_thr(m2, 0.5) #to turn it logical
caim2[!m2] <- 0

plotRGB(normalize_minmax(caim2)*255)

## End(Not run)

Description

Reprojects a fisheye (hemispherical) image into a panoramic view using a cylindrical projection. The output is standardized so that rows correspond to zenith angle bands and columns to azimuthal sectors.

Usage

fisheye_to_pano(r, z, a, fun = mean, angle_width = 1)

Arguments

Details

This function computes a cylindrical projection by aggregating pixel values according to their zenith and azimuth angles. Internally, it creates a segmentation grid with skygrid_segmentation() and applies extract_feature() to compute a summary statistic (e.g., mean) of pixel values within each cell.

Value

terra::SpatRaster with rows representing zenith angle bands and columns representing azimuthal sectors. The number of layers and names matches that of the input r.

Note

An early version of this function was used in Díaz et al. (2021).

References

Díaz GM, Negri PA, Lencinas JD (2021). “Toward making canopy hemispherical photography independent of illumination conditions: A deep-learning-based approach.” Agricultural and Forest Meteorology, 296, 108234. doi:10.1016/j.agrformet.2020.108234.

Examples

## Not run: 
caim <- read_caim()
z <- zenith_image(ncol(caim), lens())
a <- azimuth_image(z)
pano <- fisheye_to_pano(caim, z, a)
plotRGB(pano %>% normalize_minmax() %>% multiply_by(255))

## End(Not run)

Description

Fit the CIE sky model to data sampled from a canopy photograph using general-purpose optimization.

Usage

fit_cie_model(
  rr,
  sun_angles,
  custom_sky_coef = NULL,
  std_sky_no = NULL,
  general_sky_type = NULL,
  optim_zenith_dn = FALSE,
  method = c("Nelder-Mead", "BFGS", "CG", "SANN")
)

Arguments

Details

The method is based on Lang et al. (2010). For best results, the input data should show a linear relation between digital numbers and the amount of light reaching the sensor. See extract_radiometry() and read_caim_raw() for details. As a compromise solution, invert_gamma_correction() can be used.

Value

List with the following components:

rr

The input rr with an added pred column in sky_points, containing predicted values.

opt_result

List returned by stats::optim().

coef

Numeric vector of length five. CIE model coefficients.

sun_angles

Numeric vector of length two. Sun zenith and azimuth (degrees).

method

Character vector of length one. Optimization method used.

start

Numeric vector of length five. Starting parameters.

metric

Numeric value. Mean squared deviation as in Gauch et al. (2003).

Digitizing sky points with ImageJ

The point selection tool of ‘ImageJ’ software can be used to manually digitize points and create a CSV file from which to read coordinates (see Examples). After digitizing the points on the image, this is a recommended workflow: 1. Use the dropdown menu Analyze > Measure to open the Results window. 2. Use File > Save As... to obtain the CSV file.

Use this code to create the input sky_points from ImageJ data:

sky_points <- read.csv(path)
sky_points <- sky_points[c("Y", "X")]
colnames(sky_points) <- c("row", "col")

Digitizing sky points with QGIS

To use the QGIS software to manually digitize points, drag and drop the image in an empty project, create an new vector layer, digitize points manually, save the editions, and close the project.

To create the new vector layer, this is a recommended workflow:

1. Fo to the dropdown menu Layer > Create Layer > New Geopackage Layer...

2. Choose "point" in the Geometry type dropdown list

3. Make sure the CRS is EPSG:7589.

4. Click on the Toogle Editing icon

5. Click on the Add Points Feature icon.

Use this code to create the input sky_points from QGIS data:

sky_points <- terra::vect(path)
sky_points <- terra::extract(caim, sky_points, cells = TRUE)
sky_points <- terra::rowColFromCell(caim, sky_points$cell) %>% as.data.frame()
colnames(sky_points) <- c("row", "col")

References

Gauch HG, Hwang JTG, Fick GW (2003). “Model evaluation by comparison of model‐based predictions and measured values.” Agronomy Journal, 95(6), 1442–1446. ISSN 1435-0645. doi:10.2134/agronj2003.1442.

Lang M, Kuusk A, Mõttus M, Rautiainen M, Nilson T (2010). “Canopy gap fraction estimation from digital hemispherical images using sky radiance models and a linear conversion method.” Agricultural and Forest Meteorology, 150(1), 20–29. doi:10.1016/j.agrformet.2009.08.001.

Examples

## Not run: 
caim <- read_caim()
z <- zenith_image(ncol(caim), lens())
a <- azimuth_image(z)

# Manual method following Lang et al. (2010)
path <- system.file("external/sky_points.csv",
                    package = "rcaiman")
sky_points <- read.csv(path)
sky_points <- sky_points[c("Y", "X")]
colnames(sky_points) <- c("row", "col")
head(sky_points)
plot(caim$Blue)
points(sky_points$col, nrow(caim) - sky_points$row, col = 2, pch = 10)

# x11()
# plot(caim$Blue)
# sun_angles <- click(c(z, a), 1) %>% as.numeric()
sun_angles <- c(z = 49.5, a = 27.4) #taken with above lines then hardcoded

sun_row_col <- row_col_from_zenith_azimuth(z, a,
                                           sun_angles["z"],
                                           sun_angles["a"])
points(sun_row_col[2], nrow(caim) - sun_row_col[1],
       col = "yellow", pch = 8, cex = 3)

rr <- extract_rr(caim$Blue, z, a, sky_points)

set.seed(7)
model <- fit_cie_model(rr, sun_angles,
                           general_sky_type = "Clear")

plot(model$rr$sky_points$rr, model$rr$sky_points$pred)
abline(0,1)
lm(model$rr$sky_points$pred~model$rr$sky_points$rr) %>% summary()

sky <- cie_image(z, a, model$sun_angles, model$coef) * model$rr$zenith_dn

plot(sky)
ratio <- caim$Blue/sky
plot(ratio)
plot(ratio > 1.05)
plot(ratio > 1.15)

## End(Not run)

Description

Fit a polynomial model to predict relative radiance from spherical coordinates using data sampled from a canopy photograph.

Usage

fit_coneshaped_model(sky_points, method = "zenith_n_azimuth")

Arguments

Details

This model requires only sky_points, making it useful in workflows where sun position cannot be reliably estimated, such as in apply_by_direction(). Otherwise, fit_cie_model() is a better choice.

Depending on method, it can fit:

A zenith-only quadratic model

$sDN = a + b\theta + c\theta^2$

A zenith-plus-azimuth model, adding sinusoidal terms

$sDN = a + b\theta + c\theta^2 + d\sin(\phi) + e\cos(\phi)$

See Díaz and Lencinas (2018) for details on the full model.

Value

List with the following components:

fun

Function taking zenith and azimuth (degrees) and returning predicted relative radiance.

model

lm object fitted by stats::lm().

Returns NULL (with a warning) if the number of input points is fewer than 20.

References

Díaz GM, Lencinas JD (2018). “Model-based local thresholding for canopy hemispherical photography.” Canadian Journal of Forest Research, 48(10), 1204–1216. doi:10.1139/cjfr-2018-0006.

Examples

## Not run: 
caim <- read_caim()
z <- zenith_image(ncol(caim), lens())
a <- azimuth_image(z)
m <- !is.na(z)
r <- caim$Blue

bin <- binarize_by_region(r, ring_segmentation(z, 15), "thr_isodata") &
  select_sky_region(z, 0, 88)

g <- skygrid_segmentation(z, a, 10, first_ring_different = TRUE)
sky_points <- sample_sky_points(r, bin, g, dist_to_black = 3)
plot(bin)
points(sky_points$col, nrow(caim) - sky_points$row, col = 2, pch = 10)
rr <- extract_rr(r, z, a, sky_points)

model <- fit_coneshaped_model(rr$sky_points)
summary(model$model)
sky_cs <- model$fun(z, a) * rr$zenith_dn
plot(sky_cs)

z_mini <- zenith_image(50, lens())
sky_cs <- model$fun(z_mini, azimuth_image(z_mini))
persp(sky_cs, theta = 90, phi = 20)


## End(Not run)

Description

Fits a trend surface to sky digital numbers using spatial::surf.ls() as the computational workhorse.

Usage

fit_trend_surface(sky_points, r, np = 6, col_id = "dn", extrapolate = FALSE)

Arguments

Details

This function models the variation in digital numbers across the sky dome by fitting a polynomial surface in Cartesian space. It is intended to capture smooth large-scale gradients and is more effective when called via apply_by_direction().

Value

List with named elements:

raster

terra::SpatRaster containing the fitted surface.

model

object of class trls returned by spatial::surf.ls().

r2

numeric value giving the coefficient of determination (R$^2$) of the fit.

References

There are no references for Rd macro ⁠\insertAllCites⁠ on this help page.

Examples

## Not run: 
caim <- read_caim()
z <- zenith_image(ncol(caim), lens())
a <- azimuth_image(z)
m <- !is.na(z)
r <- caim$Blue

bin <- binarize_by_region(r, ring_segmentation(z, 15), "thr_isodata") &
  select_sky_region(z, 0, 88)

g <- skygrid_segmentation(z, a, 10, first_ring_different = TRUE)
sky_points <- sample_sky_points(r, bin, g, dist_to_black = 3)
plot(bin)
points(sky_points$col, nrow(caim) - sky_points$row, col = 2, pch = 10)
sky_points <- extract_dn(r, sky_points, use_window = TRUE)

sky_s <- fit_trend_surface(sky_points, r, np = 4, col_id = 3,
                           extrapolate = TRUE)
plot(sky_s$raster)
binarize_with_thr(r/sky_s$raster, 0.5) %>% plot()

sky_s <- fit_trend_surface(sky_points, r, np = 6, col_id = 3,
                           extrapolate = FALSE)
plot(sky_s$raster)
binarize_with_thr(r/sky_s$raster, 0.5) %>% plot()

## End(Not run)

Description

Grow black pixels in a binary mask using a kernel of $3 \times 3$ pixels iteratively. Useful to reduce errors associated with inter-class borders.

Usage

grow_black(bin, dist_to_black)

Arguments

Details

Expands the regions with value FALSE (typically rendered as black) in a binary image by applying a square-shaped buffer. Any white pixels (value TRUE) within a distance equal to or less than dist_to_black from a black pixel will be turned black.

Value

Logical terra::SpatRaster with the same dimensions as bin. Compared to the input bin, black regions (FALSE) have been expanded by the specified buffer distance.

Examples

## Not run: 
r <- read_caim()
bin <- binarize_with_thr(r$Blue, thr_isodata(r$Blue[]))
plot(bin)
bin <- grow_black(bin, 11)
plot(bin)

## End(Not run)

Description

Read and write legacy files from HSP (HemiSPherical Project Manager) projects to interoperate with existing workflows. Intended for legacy support; not required when working fully within rcaiman.

Usage

hsp_read_manual_input(path_to_HSP_project, img_name)

hsp_read_opt_sky_coef(path_to_HSP_project, img_name)

hsp_write_sky_points(sky_points, path_to_HSP_project, img_name)

hsp_write_sun_coord(sun_row_col, path_to_HSP_project, img_name)

Arguments

Value

See Functions

About HSP software

HSP was introduced in (Lang et al. 2013) and is based on the method presented in Lang et al. (2010)). It runs exclusively on Windows. HSP stores pre-processed images as PGM files in the manipulate subfolder of each project (itself inside the projects folder).

Functions

hsp_read_manual_input()

read sky marks and sun position defined manually within an HSP project; returns a named list with components weight, max_points, angle, point_radius, sun_row_col, sky_points, zenith_dn, masks, and marked_areas.

hsp_read_opt_sky_coef()

read optimized CIE sky coefficients from an HSP project; returns a numeric vector of length five.

hsp_write_sky_points()

write a file with sky point coordinates compatible with HSP; creates a file on disk.

hsp_write_sun_coord()

write a file with solar disk coordinates compatible with HSP; creates a file on disk.

References

Lang M, Kodar A, Arumäe T (2013). “Restoration of above canopy reference hemispherical image from below canopy measurements for plant area index estimation in forests.” Forestry Studies, 59(1), 13–27. doi:10.2478/fsmu-2013-0008.

Lang M, Kuusk A, Mõttus M, Rautiainen M, Nilson T (2010). “Canopy gap fraction estimation from digital hemispherical images using sky radiance models and a linear conversion method.” Agricultural and Forest Meteorology, 150(1), 20–29. doi:10.1016/j.agrformet.2009.08.001.

Examples

## Not run: 
# NOTE: assumes the working directory is the HSP project folder (e.g., an RStudio project).

# From HSP to R in order to compare ---------------------------------------
r <- read_caim("manipulate/IMG_1013.pgm")
z <- zenith_image(ncol(r), lens())
a <- azimuth_image(z)
manual_input <- read_manual_input(".", "IMG_1013")
sun_row_col <- manual_input$sun_row_col
sun_angles <- zenith_azimuth_from_row_col(
  z, a,
  sun_row_col[1],
  sun_row_col[2]
)
sun_angles <- as.vector(sun_angles)

sky_points <- manual_input$sky_points
rr <- extract_rr(r, z, a, sky_points)
model <- fit_cie_model(rr, sun_angles)
sky <- cie_image(
  z, a,
  model$sun_angles,
  model$coef
) * model$rr$zenith_dn
plot(r / sky)

r <- read_caim("manipulate/IMG_1013.pgm")
sky_coef <- read_opt_sky_coef(".", "IMG_1013")
sky_m <- cie_image(z, a, sun_angles, sky_coef)
sky_m <- cie_sky_manual * manual_input$zenith_dn
plot(r / sky_m)

# From R to HSP  ----------------------------------------------------------
r <- read_caim("manipulate/IMG_1014.pgm")
z <- zenith_image(ncol(r), lens())
a <- azimuth_image(z)
m <- !is.na(z)
g <- skygrid_segmentation(z, a, 10)
bin <- binarize_with_thr(caim$Blue, thr_isodata(caim$Blue[m]))
bin <- select_sky_region(z, 0, 85) & bin

sun_angles <- estimate_sun_angles(r, z, a, bin, g)
sun_row_col <- row_col_from_zenith_azimuth(
  z, a,
  sun_angles["z"],
  sun_angles["a"]
) %>% as.numeric()
write_sun_coord(sun_row_col, ".", "IMG_1014")

sky_points <- sample_sky_points(r, bin, g)
write_sky_points(sky_points, ".", "IMG_1014")

## End(Not run)

Description

Interpolate values from canopy photographs using inverse distance weighting (IDW) with k-nearest neighbors in image (planar) coordinates. A radius limits neighbor search.

Usage

interpolate_planar(
  sky_points,
  r,
  k = 3,
  p = 2,
  rmax = 200,
  col_id = "dn",
  rigorous = FALSE,
  engine = "lidR"
)

Arguments

Details

Defaults follow Lang et al. (2013). Note that rmax is given in pixels but intended to approximate 15–20 deg in angular terms. Therefore, this value needs fine-tuning based on image resolution and lens projection. For best results, the interpolated quantity should be linearly related to scene radiance; see extract_radiometry() and read_caim_raw(). For JPEG images, consider invert_gamma_correction() to reverse gamma encoding.

Value

Numeric terra::SpatRaster with one layer and the same geometry as r.

Note

No consistency checks are performed to ensure that sky_points and r are geometrically compatible. Incorrect combinations may lead to invalid outputs.

References

Lang M, Kodar A, Arumäe T (2013). “Restoration of above canopy reference hemispherical image from below canopy measurements for plant area index estimation in forests.” Forestry Studies, 59(1), 13–27. doi:10.2478/fsmu-2013-0008.

Examples

## Not run: 
caim <- read_caim()
z <- zenith_image(ncol(caim), lens())
a <- azimuth_image(z)
m <- !is.na(z)
r <- caim$Blue

bin <- binarize_by_region(r, ring_segmentation(z, 15), "thr_isodata") &
  select_sky_region(z, 0, 88)

g <- skygrid_segmentation(z, a, 10)
sky_points <- sample_sky_points(r, bin, g, dist_to_black = 3)
plot(bin)
points(sky_points$col, nrow(caim) - sky_points$row, col = 2, pch = 10)
sky_points <- extract_dn(r, sky_points)

sky <- interpolate_planar(sky_points, r, col_id = 3)
plot(sky)
plot(r/sky)

## End(Not run)

Description

Interpolate values from canopy photographs using inverse distance weighting (IDW) with k-nearest neighbors, computing distances in spherical coordinates that map the sky vault. Optionally blend with a model surface to fill voids.

Usage

interpolate_spherical(
  sky_points,
  z,
  a,
  filling_source = NULL,
  w = 1,
  k = 3,
  p = 2,
  angular_radius = 20,
  rule = "any",
  size = 50
)

Arguments

Details

Distances are great-circle distances on the sky vault (see calc_spherical_distance()). When filling_source is provided, local IDW estimates are blended with that surface following Eq. 6 in Lang et al. (2010). For efficiency, interpolation runs on a temporary low-resolution grid of size size.

Value

Numeric terra::SpatRaster with one layer of interpolated values and the geometry of z.

Note

This function assumes that sky_points and the terra::SpatRaster inputs are spatially aligned and share the same geometry. No checks are performed to enforce this.

References

Lang M, Kuusk A, Mõttus M, Rautiainen M, Nilson T (2010). “Canopy gap fraction estimation from digital hemispherical images using sky radiance models and a linear conversion method.” Agricultural and Forest Meteorology, 150(1), 20–29. doi:10.1016/j.agrformet.2009.08.001.

Examples

## Not run: 
caim <- read_caim()
z <- zenith_image(ncol(caim), lens())
a <- azimuth_image(z)

# Manual method following Lang et al. (2010)
path <- system.file("external/sky_points.csv",
                    package = "rcaiman")
sky_points <- read.csv(path)
sky_points <- sky_points[c("Y", "X")]
colnames(sky_points) <- c("row", "col")
head(sky_points)
plot(caim$Blue)
points(sky_points$col, nrow(caim) - sky_points$row, col = 2, pch = 10)

# x11()
# plot(caim$Blue)
# sun_angles <- click(c(z, a), 1) %>% as.numeric()
sun_angles <- c(z = 49.5, a = 27.4) #taken with above lines then hardcoded

sun_row_col <- row_col_from_zenith_azimuth(z, a,
                                           sun_angles["z"],
                                           sun_angles["a"])
points(sun_row_col[2], nrow(caim) - sun_row_col[1],
       col = "yellow", pch = 8, cex = 3)

rr <- extract_rr(caim$Blue, z, a, sky_points)

set.seed(7)
model <- fit_cie_model(rr, sun_angles,
                       general_sky_type = "Clear")

sky_cie <- cie_image(z, a, model$sun_angles, model$coef)

sky_rr <- interpolate_spherical(rr$sky_points, z, a,
                                filling_source = sky_cie,
                                w = 1,
                                k = 10,
                                p = 2,
                                angular_radius = 20,
                                rule = "any",
                                size = 50)
plot(caim$Blue/sky_rr/rr$zenith_dn)

## End(Not run)

Description

Approximates the inversion of the gamma encoding applied to JPEG images.

Usage

invert_gamma_correction(dn, gamma = 2.2)

Arguments

Details

Digital cameras typically encode images using the sRGB color space, which applies a non-linear transformation—commonly referred to as gamma correction—to the sensor's linear luminance response. This function applies a power transformation to approximate the inverse of that encoding, restoring a response closer to linear.

Value

Same properties as dn, with values adjusted by inverse gamma correction and rescaled to the range $[0, 1]$.

References

There are no references for Rd macro ⁠\insertAllCites⁠ on this help page.

Examples

## Not run: 
path <- system.file("external/APC_0836.jpg", package = "rcaiman")
caim <- read_caim(path)
z <- zenith_image(2132,  lens("Olloclip"))
a <- azimuth_image(z)
zenith_colrow <- c(1063, 771)

caim <- expand_noncircular(caim, z, zenith_colrow)
m <- !is.na(caim$Red) & !is.na(z)
caim[!m] <- 0

bin <- binarize_with_thr(caim$Blue, thr_isodata(caim$Blue[m]))

display_caim(caim$Blue, bin)

caim <- invert_gamma_correction(caim, 2.2)

## End(Not run)

Description

Identify outlying values in a numeric vector using a robust MAD–median rule and return a logical mask indicating which elements are outliers.

Usage

is_outlier(x, y = NULL, laxity = 2, cutoff_side = "both")

Arguments

Details

Apply the robust criterion from Leys et al. (2013) to detect statistical outliers. The function returns a logical vector of the same length as x with TRUE for detected outliers.

An observation $x_i$ is considered outlying if it satisfies:

$M - \mathrm{laxity} \times \mathrm{MAD} < x_i < M + \mathrm{laxity} \times \mathrm{MAD}$

where $M$ is the median of x and $MAD$ is its median absolute deviation. The laxity parameter multiplies MAD to set the classification threshold.

The cutoff_side argument controls which tail(s) are evaluated: "both" (default) applies a two-sided test, "left" applies the left-tail test only, and "right" applies the right-tail test only.

If the MAD is zero (no dispersion), all values are flagged as non-outliers. NA values in x are treated as non-outliers.

Value

Logical vector of the same length as x. TRUE indicates the outliers.

References

Leys C, Ley C, Klein O, Bernard P, Licata L (2013). “Detecting outliers: Do not use standard deviation around the mean, use absolute deviation around the median.” Journal of Experimental Social Psychology, 49(4), 764–766. ISSN 0022-1031. doi:10.1016/j.jesp.2013.03.013.

Examples

x <- c(rnorm(100), 10, 12)
i <- is_outlier(x, x, laxity = 2.5, cutoff_side = "both")
x[i]  # values flagged as outliers
x[!i] # filter outliers

Description

Retrieve projection coefficients and field-of-view (FOV, deg) for known lenses. Coefficients expect zenith angle in radians and return relative radius.

Usage

lens(type = "equidistant", return = "coef")

Arguments

Details

In upward-looking leveled hemispherical photography, the zenith corresponds to the center of a circular image whose perimeter represents the horizon, assuming a lens with a 180° field of view. The relative radius is the radial distance to a given point, expressed as a fraction of the maximum radius (i.e., the horizon). The equidistant projection model, also called polar projection, is the standard reference model, characterized by a linear relationship between zenith angle and relative radius.

Real lenses deviate from ideal projections. Following Hemisfer software, this package models deviations with polynomial functions. All angular values are in radians.

Currently available lenses:

"equidistant"

Ideal equidistant projection (Schneider et al. 2009).

"Nikkor_10.5mm"

AF DX Fisheye Nikkor 10.5mm f/2.8G ED (Pekin and Macfarlane 2009).

"Nikon_FCE9"

Nikon FC-E9 converter (Díaz et al. 2024).

"Olloclip"

Auxiliary lens for mobile devices manufactured by Olloclip (Díaz et al. 2024).

"Nikkor_8mm"

AF–S Fisheye Nikkor 8–15mm f/3.5–4.5E ED (Díaz et al. 2024).

See calibrate_lens() for fitting new projection functions.

Value

Numeric vector. Depends on return:

"coef"

Polynomial coefficients of the projection function (relative radius as a function of theta, radians).

"fov"

numeric vector of length one. Maximum field of view (deg).

References

Díaz GM, Lang M, Kaha M (2024). “Simple calibration of fisheye lenses for hemipherical photography of the forest canopy.” Agricultural and Forest Meteorology, 352, 110020. ISSN 0168-1923. doi:10.1016/j.agrformet.2024.110020.

Pekin B, Macfarlane C (2009). “Measurement of crown cover and leaf area index using digital cover photography and its application to remote sensing.” Remote Sensing, 1(4), 1298–1320. doi:10.3390/rs1041298.

Schneider D, Schwalbe E, Maas H (2009). “Validation of geometric models for fisheye lenses.” ISPRS Journal of Photogrammetry and Remote Sensing, 64(3), 259–266. doi:10.1016/j.isprsjprs.2009.01.001.

Examples

lens("Nikon_FCE9")
lens("Nikon_FCE9", return = "fov")

.fp <- function(theta, lens_coef) {
  x <- lens_coef[1:5]
  x[is.na(x)] <- 0
  for (i in 1:5) assign(letters[i], x[i])
  a * theta + b * theta^2 + c * theta^3 + d * theta^4 + e * theta^5
}

theta <- seq(0, pi/2, pi/180)
plot(theta, .fp(theta, lens()), type = "l", lty = 2,
      ylab = "relative radius")
lines(theta, .fp(theta, lens("Nikon_FCE9")))

Description

Rescale numeric or raster data from an expected range to the range $[0, 1]$.

Usage

normalize_minmax(r, mn = NULL, mx = NULL, clip = FALSE)

Arguments

Value

Same properties as r, with values rescaled to the range $[0, 1]$ if mn and mx match the range of r or extend beyond it. If clip = TRUE, values will be within $[0, 1]$ even if this implies data loss.