Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
194 changes: 97 additions & 97 deletions src/compute.rs
Original file line number Diff line number Diff line change
Expand Up @@ -116,15 +116,19 @@ impl Warper {
*v = result;
FoldWhile::Continue(Ok(()))
} else {
FoldWhile::Done(Err(WarperError::NanError))
FoldWhile::Done(Err(WarperError::WarpingError))
}
})
.into_inner()?;

Ok(target_raster)
}

/// This variant throws an error if there's any NaN in the data or result is not finite
/// This variant throws an error if there's any NaN in the (used) input data or result is not finite
///
/// Note: the input NaN check checks only values that are actually used in computations,
/// so you can pass an input raster with NaNs and get `Ok` result as long as those values
/// do not appear in any computation kernel.
pub fn warp_reject_nodata<'a, A: Into<ArrayView2<'a, f64>>>(
&self,
source_raster: A,
Expand Down Expand Up @@ -154,7 +158,7 @@ impl Warper {
let value = values[[j, i]];

if value.is_nan() {
return FoldWhile::Done(Err(WarperError::WarpingError));
return FoldWhile::Done(Err(WarperError::NanError));
}
let x_weight = intr.x_weights[i];
inner_weight_accum += x_weight;
Expand Down Expand Up @@ -283,13 +287,15 @@ impl Warper {
}

/// This implementation catches warp operation producing NaNs (that is nans resulting from computation error
/// not those resulting from nodata), but does not early return.
/// not those resulting from nodata)
#[cfg(feature = "multithreading")]
#[cfg_attr(docsrs, doc(cfg(feature = "multithreading")))]
pub fn warp_ignore_nodata_parallel<'a, A: Into<ArrayView2<'a, f64>>>(
&self,
source_raster: A,
) -> Result<Array2<f64>, WarperError> {
use ndarray::parallel::prelude::{IntoParallelIterator, ParallelIterator};

let source_raster: ArrayView2<f64> = source_raster.into();

self.validate_source_raster_shape(&source_raster)?;
Expand All @@ -298,53 +304,50 @@ impl Warper {

Zip::from(&mut target_raster)
.and(&self.internals)
.par_fold(
|| Ok(()),
|_, v, intr| {
let values = source_raster.slice(s![
(intr.anchor_idx.1 - 1)..(intr.anchor_idx.1 + 3),
(intr.anchor_idx.0 - 1)..(intr.anchor_idx.0 + 3)
]);

let mut weight_accum = 0.0;
let mut result_accum = 0.0;

for j in 0..4 {
let mut inner_weight_accum = 0.0;
let mut inner_result_accum = 0.0;

for i in 0..4 {
let value = values[[j, i]];

if !value.is_nan() {
let x_weight = intr.x_weights[i];
inner_weight_accum += x_weight;
inner_result_accum += x_weight * value;
}
}
.into_par_iter()
.try_for_each(|(v, intr)| {
let values = source_raster.slice(s![
(intr.anchor_idx.1 - 1)..(intr.anchor_idx.1 + 3),
(intr.anchor_idx.0 - 1)..(intr.anchor_idx.0 + 3)
]);

let mut weight_accum = 0.0;
let mut result_accum = 0.0;

let y_weight = intr.y_weights[j];
for j in 0..4 {
let mut inner_weight_accum = 0.0;
let mut inner_result_accum = 0.0;

weight_accum += inner_weight_accum * y_weight;
result_accum += inner_result_accum * y_weight;
}
for i in 0..4 {
let value = values[[j, i]];

if (weight_accum).abs() < f64::EPSILON {
*v = f64::NAN;
return Ok(());
if !value.is_nan() {
let x_weight = intr.x_weights[i];
inner_weight_accum += x_weight;
inner_result_accum += x_weight * value;
}
}

let result = result_accum / weight_accum;
let y_weight = intr.y_weights[j];

if result.is_finite() {
*v = result;
Ok(())
} else {
Err(WarperError::WarpingError)
}
},
std::result::Result::and,
)?;
weight_accum += inner_weight_accum * y_weight;
result_accum += inner_result_accum * y_weight;
}

if (weight_accum).abs() < f64::EPSILON {
*v = f64::NAN;
return Ok(());
}

let result = result_accum / weight_accum;

if result.is_finite() {
*v = result;
Ok(())
} else {
Err(WarperError::WarpingError)
}
})?;

Ok(target_raster)
}
Expand All @@ -356,47 +359,47 @@ impl Warper {
&self,
source_raster: A,
) -> Result<Array2<f64>, WarperError> {
use ndarray::parallel::prelude::{IntoParallelIterator, ParallelIterator};

let source_raster: ArrayView2<f64> = source_raster.into();

self.validate_source_raster_shape(&source_raster)?;

Zip::from(&source_raster).par_fold(
|| Ok(()),
|_, v| {
Zip::from(&source_raster)
.into_par_iter()
.try_for_each(|(v,)| {
if v.is_nan() {
Err(WarperError::NanError)
} else {
Ok(())
}
},
std::result::Result::and,
)?;
})?;

let target_raster = self.warp_unchecked_parallel(source_raster);

Zip::from(&target_raster).par_fold(
|| Ok(()),
|_, v| {
Zip::from(&target_raster)
.into_par_iter()
.try_for_each(|(v,)| {
if v.is_finite() {
Ok(())
} else {
Err(WarperError::WarpingError)
}
},
std::result::Result::and,
)?;
})?;

Ok(target_raster)
}

/// This implementation catches warp operation producing NaNs (that is nans resulting from computation error
/// not those resulting from nodata), but does not early return.
/// not those resulting from nodata)
#[cfg(feature = "multithreading")]
#[cfg_attr(docsrs, doc(cfg(feature = "multithreading")))]
pub fn warp_discard_nodata_parallel<'a, A: Into<ArrayView2<'a, f64>>>(
&self,
source_raster: A,
) -> Result<Array2<f64>, WarperError> {
use ndarray::parallel::prelude::{IntoParallelIterator, ParallelIterator};

let source_raster: ArrayView2<f64> = source_raster.into();

self.validate_source_raster_shape(&source_raster)?;
Expand All @@ -405,50 +408,47 @@ impl Warper {

Zip::from(&mut target_raster)
.and(&self.internals)
.par_fold(
|| Ok(()),
|_, v, intr| {
let values = source_raster.slice(s![
(intr.anchor_idx.1 - 1)..(intr.anchor_idx.1 + 3),
(intr.anchor_idx.0 - 1)..(intr.anchor_idx.0 + 3)
]);

let mut weight_accum = 0.0;
let mut result_accum = 0.0;

for j in 0..4 {
let mut inner_weight_accum = 0.0;
let mut inner_result_accum = 0.0;

for i in 0..4 {
let value = values[[j, i]];

if value.is_nan() {
*v = f64::NAN;
return Ok(());
}
let x_weight = intr.x_weights[i];
inner_weight_accum += x_weight;
inner_result_accum += x_weight * value;
}
.into_par_iter()
.try_for_each(|(v, intr)| {
let values = source_raster.slice(s![
(intr.anchor_idx.1 - 1)..(intr.anchor_idx.1 + 3),
(intr.anchor_idx.0 - 1)..(intr.anchor_idx.0 + 3)
]);

let y_weight = intr.y_weights[j];
let mut weight_accum = 0.0;
let mut result_accum = 0.0;

weight_accum += inner_weight_accum * y_weight;
result_accum += inner_result_accum * y_weight;
}
for j in 0..4 {
let mut inner_weight_accum = 0.0;
let mut inner_result_accum = 0.0;

let result = result_accum / weight_accum;
for i in 0..4 {
let value = values[[j, i]];

if result.is_finite() {
*v = result;
Ok(())
} else {
Err(WarperError::WarpingError)
if value.is_nan() {
*v = f64::NAN;
return Ok(());
}
let x_weight = intr.x_weights[i];
inner_weight_accum += x_weight;
inner_result_accum += x_weight * value;
}
},
std::result::Result::and,
)?;

let y_weight = intr.y_weights[j];

weight_accum += inner_weight_accum * y_weight;
result_accum += inner_result_accum * y_weight;
}

let result = result_accum / weight_accum;

if result.is_finite() {
*v = result;
Ok(())
} else {
Err(WarperError::WarpingError)
}
})?;

Ok(target_raster)
}
Expand Down
17 changes: 10 additions & 7 deletions src/precompute.rs
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,10 @@ pub fn precompute_ixs_jys_parallel<SP: Projection, TP: Projection>(
source_bounds: &RasterBounds<SP, SourceXYPair>,
target_bounds: &RasterBounds<TP, TargetXYPair>,
) -> Result<Array2<IXJYPair>, WarperError> {
use ndarray::Zip;
use ndarray::{
Zip,
parallel::prelude::{IntoParallelIterator, ParallelIterator},
};

let tgt_ul_edge_corner = SourceXYPair {
x: 0.5f64.mul_add(-target_bounds.spacing.x, target_bounds.min.x),
Expand Down Expand Up @@ -98,17 +101,17 @@ pub fn precompute_ixs_jys_parallel<SP: Projection, TP: Projection>(
)
});

Zip::from(&precomputed_coords).par_fold(
|| Ok(()),
|_, v| {
// ndarray uses into_par_iter() under the hood so we are not loosing performance
// by going to rayon here, possibly even gaining something
Zip::from(&precomputed_coords)
.into_par_iter()
.try_for_each(|(v,)| {
if !v.ix.is_finite() || !v.jy.is_finite() {
Err(WarperError::ConversionError)
} else {
Ok(())
}
},
std::result::Result::and,
)?;
})?;

Ok(precomputed_coords)
}
Expand Down
Loading