diff --git a/Cargo.toml b/Cargo.toml
index 16bb8de..b5bb80c 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -1,6 +1,6 @@
 [package]
 name = "imgproc-rs"
-version = "0.2.1"
+version = "0.2.2"
 edition = "2018"
 license = "MIT"
 description = "Image processing library for Rust"
diff --git a/README.md b/README.md
index 14b2b96..d39194c 100644
--- a/README.md
+++ b/README.md
@@ -128,7 +128,7 @@ To enable multithreading, include the `parallel` feature in your `Cargo.toml`:
 
 ```toml
 [dependencies.imgproc-rs]
-version = "0.2.1"
+version = "0.2.2"
 default-features = false
 features = ["parallel"]
 ```
diff --git a/src/colorspace.rs b/src/colorspace.rs
index 334e088..fc6bdd0 100644
--- a/src/colorspace.rs
+++ b/src/colorspace.rs
@@ -1,12 +1,12 @@
 //! A module for image colorspace conversion operations
 
-use crate::{util, math};
-use crate::util::constants::{GAMMA, SRGB_TO_XYZ_MAT, XYZ_TO_SRGB_MAT};
-use crate::image::Image;
-use crate::enums::White;
-
 use std::cmp;
 
+use crate::enums::White;
+use crate::image::Image;
+use crate::util;
+use crate::util::constants::{GAMMA, SRGB_TO_XYZ_MAT, XYZ_TO_SRGB_MAT};
+
 /// Converts an image from RGB to Grayscale
 pub fn rgb_to_grayscale(input: &Image<u8>) -> Image<u8> {
     input.map_pixels_if_alpha(|channels| {
@@ -32,11 +32,12 @@ pub fn rgb_to_grayscale_f64(input: &Image<f64>) -> Image<f64> {
 }
 
 /// Linearizes an sRGB image
-// Input: sRGB range [0, 255]
-// Output: sRGB range [0, 1] linearized
+///
+/// * Input: sRGB image with channels in range [0, 255]
+/// * Output: linearized sRGB image with channels in range [0, 1]
 pub fn linearize_srgb(input: &Image<u8>) -> Image<f64> {
     let mut lookup_table: [f64; 256] = [0.0; 256];
-    util::create_lookup_table(&mut lookup_table, |i| {
+    util::generate_lookup_table(&mut lookup_table, |i| {
         let val = i as f64;
         if val <= 10.0 {
             val / 3294.0
@@ -49,8 +50,9 @@ pub fn linearize_srgb(input: &Image<u8>) -> Image<f64> {
 }
 
 /// "Unlinearizes" a previously linearized sRGB image
-// Input: sRGB range [0, 1] linearized
-// Output: sRGB range [0, 255]
+///
+/// * Input: linearized sRGB image with channels in range [0, 1]
+/// * Output: sRGB image with channels in range [0, 255]
 pub fn unlinearize_srgb(input: &Image<f64>) -> Image<u8> {
     input.map_channels_if_alpha(|num| {
         if num <= 0.0031308 {
@@ -62,28 +64,31 @@ pub fn unlinearize_srgb(input: &Image<f64>) -> Image<u8> {
 }
 
 /// Converts an image from linearized sRGB to CIE XYZ
-// Input: sRGB range [0, 1] linearized
-// Output: CIE XYZ range [0, 1]
+///
+/// * Input: linearized sRGB image with channels in range [0, 1]
+/// * Output: CIE XYZ image with channels in range [0, 1]
 pub fn srgb_lin_to_xyz(input: &Image<f64>) -> Image<f64> {
     input.map_pixels_if_alpha(|channels| {
-        math::vector_mul(&SRGB_TO_XYZ_MAT, channels).unwrap()
+        util::vector_mul(&SRGB_TO_XYZ_MAT, channels).unwrap()
     }, |a| a)
 }
 
 /// Converts an image from CIE XYZ to linearized sRGB
-// Input: CIE XYZ range [0, 1]
-// Output: sRGB range [0, 1] linearized
+///
+/// * Input: CIE XYZ image with channels in range [0, 1]
+/// * Output: linearized sRGB image with channels in range [0, 1]
 pub fn xyz_to_srgb_lin(input: &Image<f64>) -> Image<f64> {
     input.map_pixels_if_alpha(|channels| {
-        math::vector_mul(&XYZ_TO_SRGB_MAT, channels).unwrap()
+        util::vector_mul(&XYZ_TO_SRGB_MAT, channels).unwrap()
     }, |a| a)
 }
 
 /// Converts an image from CIE XYZ to CIELAB
-// Input: CIEXYZ range [0, 1]
-// Output: CIELAB with L* channel range [0, 100] and a*, b* channels range [-128,127]
+///
+/// * Input: CIE XYZ image with channels in range [0, 1]
+/// * Output: CIELAB image with L* channel range [0, 100] and a*, b* channels range [-128, 127]
 pub fn xyz_to_lab(input: &Image<f64>, ref_white: &White) -> Image<f64> {
-    let (x_n, y_n, z_n) = util::generate_xyz_tristimulus_vals(ref_white);
+    let (x_n, y_n, z_n) = util::xyz_tristimulus_vals(ref_white);
 
     input.map_pixels_if_alpha(|channels| {
         let x = util::xyz_to_lab_fn(channels[0] * 100.0 / x_n);
@@ -97,10 +102,11 @@ pub fn xyz_to_lab(input: &Image<f64>, ref_white: &White) -> Image<f64> {
 }
 
 /// Converts an image from CIELAB to CIE XYZ
-// Input: CIELAB with L* channel range [0, 100] and a*, b* channels range [-128,127]
-// Output: CIEXYZ range [0, 1]
+///
+/// * Input: CIELAB image with L* channel range [0, 100] and a*, b* channels range [-128, 127]
+/// * Output: CIE XYZ image with channels in range [0, 1]
 pub fn lab_to_xyz(input: &Image<f64>, ref_white: &White) -> Image<f64> {
-    let (x_n, y_n, z_n) = util::generate_xyz_tristimulus_vals(ref_white);
+    let (x_n, y_n, z_n) = util::xyz_tristimulus_vals(ref_white);
 
     input.map_pixels_if_alpha(|channels| {
         let n = (channels[0] + 16.0) / 116.0;
@@ -112,8 +118,9 @@ pub fn lab_to_xyz(input: &Image<f64>, ref_white: &White) -> Image<f64> {
 }
 
 /// Converts an image from RGB to HSV
-// Input: RGB range [0, 255]
-// Output: HSV range [0, 1]
+///
+/// * Input: RGB image with channels in range [0, 255]
+/// * Output: HSV image with channels in range [0, 1]
 pub fn rgb_to_hsv(input: &Image<u8>) -> Image<f64> {
     input.map_pixels_if_alpha(|channels| {
         let max: u8 = cmp::max(cmp::max(channels[0], channels[1]), channels[2]);
@@ -150,8 +157,9 @@ pub fn rgb_to_hsv(input: &Image<u8>) -> Image<f64> {
 }
 
 /// Converts an image from HSV to RGB
-// Input: HSV range [0, 1]
-// Output: RGB range [0, 255]
+///
+/// * Input: HSV image with channels in range [0, 1]
+/// * Output: RGB image with channels in range [0, 255]
 pub fn hsv_to_rgb(input: &Image<f64>) -> Image<u8> {
     input.map_pixels_if_alpha(|channels| {
         if channels[1] == 0.0 {
@@ -178,32 +186,36 @@ pub fn hsv_to_rgb(input: &Image<f64>) -> Image<u8> {
 }
 
 /// Converts an image from sRGB to CIE XYZ
-// Input: sRGB range [0, 255] unlinearized
-// Output: CIEXYZ range [0, 1]
+///
+/// * Input: sRGB image with channels in range [0, 255]
+/// * Output: CIE XYZ image with channels in range [0, 1]
 pub fn srgb_to_xyz(input: &Image<u8>) -> Image<f64> {
     let linearized = linearize_srgb(input);
     srgb_lin_to_xyz(&linearized)
 }
 
 /// Converts an image from CIE XYZ to sRGB
-// Input: CIEXYZ range [0, 1]
-// Output: sRGB range [0, 255] unlinearized
+///
+/// * Input: CIE XYZ image with channels in range [0, 1]
+/// * Output: sRGB image with channels in range [0, 255]
 pub fn xyz_to_srgb(input: &Image<f64>) -> Image<u8> {
     let srgb = xyz_to_srgb_lin(input);
     unlinearize_srgb(&srgb)
 }
 
 /// Converts an image from sRGB to CIELAB
-// Input: sRGB range [0, 255] unlinearized
-// Output: CIELAB with L* channel range [0, 100] and a*, b* channels range [-128,127]
+///
+/// * Input: sRGB image with channels in range [0, 255]
+/// * Output: CIELAB image with L* channel range [0, 100] and a*, b* channels range [-128, 127]
 pub fn srgb_to_lab(input: &Image<u8>, ref_white: &White) -> Image<f64> {
     let xyz = srgb_to_xyz(input);
     xyz_to_lab(&xyz, ref_white)
 }
 
 /// Converts an image from CIELAB to sRGB
-// Input: CIELAB with L* channel range [0, 100] and a*, b* channels range [-128,127]
-// Output: sRGB range [0, 255] unlinearized
+///
+/// * Input: CIELAB image with L* channel range [0, 100] and a*, b* channels range [-128,127]
+/// * Output: sRGB image with channels in range [0, 255]
 pub fn lab_to_srgb(input: &Image<f64>, ref_white: &White) -> Image<u8> {
     let xyz = lab_to_xyz(input, ref_white);
     xyz_to_srgb(&xyz)
diff --git a/src/convert.rs b/src/convert.rs
index f1afbfb..c420ffb 100644
--- a/src/convert.rs
+++ b/src/convert.rs
@@ -2,14 +2,12 @@
 
 use crate::image::Image;
 use crate::error::ImgProcResult;
-use crate::error;
 
 /// Scales channels from range 0.0 to `current_max` to range 0.0 to `scaled_max`
-pub fn scale_channels(input: &Image<f64>, current_max: f64, scaled_max: f64) -> ImgProcResult<Image<f64>> {
-    error::check_non_neg(current_max, "current_max")?;
-    error::check_non_neg(scaled_max, "scaled_max")?;
-
-    Ok(input.map_channels(|channel| (channel / current_max * scaled_max)))
+pub fn scale_channels(input: &Image<f64>, current_min: f64, scaled_min: f64, current_max: f64, scaled_max: f64) -> ImgProcResult<Image<f64>> {
+    Ok(input.map_channels(|channel| {
+        (channel - current_min) / (current_max - current_min) * (scaled_max - scaled_min) + scaled_min
+    }))
 }
 
 /// Converts an `Image<f64>` with channels in range 0 to `scale` to an `Image<u8>` with channels
diff --git a/src/error/messages.rs b/src/error/messages.rs
index 0511f6b..e6227c3 100644
--- a/src/error/messages.rs
+++ b/src/error/messages.rs
@@ -1,5 +1,5 @@
 use crate::error::{ImgProcResult, ImgProcError};
-use crate::image::Number;
+use crate::image::{Number, Image, BaseImage};
 
 pub(crate) fn check_channels(channels: u8, len: usize) {
     if channels != len as u8 {
@@ -57,8 +57,16 @@ pub(crate) fn check_square(val: f64, name: &str) -> ImgProcResult<()> {
     Ok(())
 }
 
-pub(crate) fn check_grayscale(channels: u8, alpha: bool) -> ImgProcResult<()> {
-    if (alpha && channels != 2) || (!alpha && channels != 1) {
+pub(crate) fn check_in_range<T: Number>(val: T, min: T, max: T, name: &str) -> ImgProcResult<()> {
+    if val < min || val > max {
+        return Err(ImgProcError::InvalidArgError(format!("{} must be between {} and {} (inclusive)", name, min, max)));
+    }
+
+    Ok(())
+}
+
+pub(crate) fn check_grayscale<T: Number>(input: &Image<T>) -> ImgProcResult<()> {
+    if (input.info().alpha && input.info().channels != 2) || (!input.info().alpha && input.info().channels != 1) {
         return Err(ImgProcError::InvalidArgError("input is not a grayscale image".to_string()));
     }
 
diff --git a/src/filter/bilateral.rs b/src/filter/bilateral.rs
index cd2f8a5..fc5c768 100644
--- a/src/filter/bilateral.rs
+++ b/src/filter/bilateral.rs
@@ -1,11 +1,11 @@
-use crate::{error, util, colorspace, math};
-use crate::enums::{Bilateral, White};
-use crate::image::{Image, BaseImage};
-use crate::error::ImgProcResult;
-
 #[cfg(feature = "rayon")]
 use rayon::prelude::*;
 
+use crate::{colorspace, error, util};
+use crate::enums::{Bilateral, White};
+use crate::error::ImgProcResult;
+use crate::image::{BaseImage, Image};
+
 /// Applies a bilateral filter using CIE LAB
 #[cfg(not(feature = "rayon"))]
 pub fn bilateral_filter(input: &Image<u8>, range: f64, spatial: f64, algorithm: Bilateral)
@@ -73,7 +73,7 @@ fn bilateral_direct_pixel(input: &Image<f64>, range: f64, spatial_mat: &[f64], s
         let mut p_curr = 0.0;
 
         for i in 0..((size * size) as usize) {
-            let g_r = math::gaussian_fn((p_in[c] - p_n[i][c]).abs(), range).unwrap();
+            let g_r = util::gaussian_fn((p_in[c] - p_n[i][c]).abs(), range).unwrap();
             let weight = spatial_mat[i] * g_r;
 
             p_curr += weight * p_n[i][c];
diff --git a/src/filter/edge.rs b/src/filter/edge.rs
new file mode 100644
index 0000000..44050aa
--- /dev/null
+++ b/src/filter/edge.rs
@@ -0,0 +1,64 @@
+////////////////////
+// Edge detection
+////////////////////
+
+use crate::{filter, error, util, convert};
+use crate::image::{Image, BaseImage};
+use crate::error::ImgProcResult;
+use crate::util::constants::{K_PREWITT_1D_VERT, K_PREWITT_1D_HORZ, K_SOBEL_1D_VERT, K_SOBEL_1D_HORZ, K_LAPLACIAN};
+
+/// Applies a separable derivative mask to a grayscale image
+pub fn derivative_mask(input: &Image<f64>, vert_kernel: &[f64], horz_kernel: &[f64]) -> ImgProcResult<Image<f64>> {
+    error::check_grayscale(input)?;
+
+    let img_x = filter::separable_filter(&input, &vert_kernel, &horz_kernel)?;
+    let img_y = filter::separable_filter(&input, &horz_kernel, &vert_kernel)?;
+
+    let mut output = Image::blank(input.info());
+
+    for i in 0..(output.info().full_size() as usize) {
+        output.set_pixel_indexed(i, &[(img_x[i][0].powf(2.0) + img_y[i][0].powf(2.0)).sqrt()]);
+    }
+
+    Ok(output)
+}
+
+/// Applies the Prewitt operator to a grayscale image
+pub fn prewitt(input: &Image<f64>) -> ImgProcResult<Image<f64>> {
+    Ok(derivative_mask(input, &K_PREWITT_1D_VERT, &K_PREWITT_1D_HORZ)?)
+}
+
+/// Applies the Sobel operator to a grayscale image
+pub fn sobel(input: &Image<f64>) -> ImgProcResult<Image<f64>> {
+    Ok(derivative_mask(input, &K_SOBEL_1D_VERT, &K_SOBEL_1D_HORZ)?)
+}
+
+/// Applies a Sobel operator with weight `weight` to a grayscale image
+pub fn sobel_weighted(input: &Image<f64>, weight: u32) -> ImgProcResult<Image<f64>> {
+    let vert_kernel = vec![1.0, weight as f64, 1.0];
+    Ok(derivative_mask(input, &vert_kernel, &K_SOBEL_1D_HORZ)?)
+}
+
+/// Applies the Laplacian operator to a grayscale image. Output contains positive
+/// and negative values - use [`normalize_laplacian()`](fn.normalize_laplacian.html) for visualization
+pub fn laplacian(input: &Image<f64>) -> ImgProcResult<Image<f64>> {
+    Ok(filter::unseparable_filter(input, &K_LAPLACIAN)?)
+}
+
+/// Applies the Laplacian of Gaussian operator using a `size x size` kernel to a grayscale image.
+/// Output contains positive and negative values - use
+/// [`normalize_laplacian()`](fn.normalize_laplacian.html) for visualization
+pub fn laplacian_of_gaussian(input: &Image<f64>, size: u32, sigma: f64) -> ImgProcResult<Image<f64>> {
+    let kernel = util::generate_log_kernel(size, sigma)?;
+    Ok(filter::unseparable_filter(input, &kernel)?)
+}
+
+/// Normalizes the result of a Laplacian or Laplacian of Gaussian operator to the range [0, 255]
+pub fn normalize_laplacian(input: &Image<f64>) -> ImgProcResult<Image<u8>> {
+    error::check_grayscale(input)?;
+
+    let min = *input.data().iter().min_by(|x, y| x.partial_cmp(y).unwrap()).unwrap();
+    let max = *input.data().iter().max_by(|x, y| x.partial_cmp(y).unwrap()).unwrap();
+
+    Ok(convert::scale_channels(&input, min, 0.0, max, 255.0)?.into())
+}
\ No newline at end of file
diff --git a/src/filter/median.rs b/src/filter/median.rs
index 917994e..279f98c 100644
--- a/src/filter/median.rs
+++ b/src/filter/median.rs
@@ -6,8 +6,8 @@ use std::cmp::Reverse;
 
 /// Applies a median filter, where each output pixel is the median of the pixels in a
 /// `(2 * radius + 1) x (2 * radius + 1)` kernel in the input image. Based on Ben Weiss' partial
-/// histogram method, using a tier radix of 2. For a detailed description, see:
-/// http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.93.1608&rep=rep1&type=pdf
+/// histogram method, using a tier radix of 2. A detailed description can be found
+/// [here](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.93.1608&rep=rep1&type=pdf).
 pub fn median_filter(input: &Image<u8>, radius: u32) -> ImgProcResult<Image<u8>> {
     let mut n_cols = (4.0 * (radius as f64).powf(2.0 / 3.0)).floor() as usize;
     if n_cols % 2 == 0 {
@@ -23,8 +23,9 @@ pub fn median_filter(input: &Image<u8>, radius: u32) -> ImgProcResult<Image<u8>>
     Ok(output)
 }
 
-/// Applies an alpha-trimmed mean filter, where each output pixel is the alpha-trimmed mean of the
-/// pixels in a `(2 * radius + 1) x (2 * radius + 1)` kernel in the input image
+/// Applies an alpha-trimmed mean filter, where each output pixel is the mean of the
+/// pixels in a `(2 * radius + 1) x (2 * radius + 1)` kernel in the input image, with the lowest
+/// `alpha / 2` pixels and the highest `alpha / 2` pixels removed.
 pub fn alpha_trimmed_mean_filter(input: &Image<u8>, radius: u32, alpha: u32) -> ImgProcResult<Image<u8>> {
     let size = 2 * radius + 1;
     error::check_even(alpha, "alpha")?;
@@ -46,13 +47,36 @@ pub fn alpha_trimmed_mean_filter(input: &Image<u8>, radius: u32, alpha: u32) ->
     Ok(output)
 }
 
+/*
+ * The PartialHistograms struct:
+ *
+ * This struct contains the partial histograms, which is a vector of an odd number of histograms
+ * determined by n_cols. The only "complete" histogram is the central histogram (located at
+ * data[n_half]), which is the histogram of the pixel values in the kernel surrounding the
+ * central pixel in the row that is being processed. Each histogram to the left and right of the
+ * central histogram is not another "complete" histogram, but rather is a histogram representing
+ * the difference between the histogram for the pixel at that location and the central histogram.
+ * As such, the values in these "partial" histograms can (and frequently will) be negative.
+ * The "complete" histogram for each non-central pixel is then just the sum of the corresponding
+ * partial histogram and the central histogram.
+ *
+ * Algorithm overview:
+ *
+ * The basic idea of this algorithm is to process a row of n_cols pixels at once using the
+ * partial histograms to efficiently compute the complete histograms for each pixel. To process the
+ * next row, the partial histograms are updated to remove the top row of pixel values from the
+ * previous kernel and add the bottom row of pixel values from the current kernel. Each set of
+ * n_cols columns in the image is processed in this fashion, using a single set of partial
+ * histograms that are updated as the current kernel slides down the image.
+ */
 #[derive(Debug, Clone)]
 struct PartialHistograms {
     data: Vec<[i32; 256]>, // The partial histograms
-    n_cols: usize,
-    n_half: usize,
-    radius: usize,
-    size: usize,
+    n_cols: usize, // The number of partial histograms, which is always odd. This also denotes the
+                   // number of columns we can process at once
+    n_half: usize, // Half the number of partial histograms, rounded down
+    radius: usize, // The radius of the kernel we are using
+    size: usize, // The number of pixels in a kernel
 }
 
 impl PartialHistograms {
@@ -69,6 +93,7 @@ impl PartialHistograms {
         }
     }
 
+    // Add or remove a row of pixels from the histograms, as indicated by the add parameter
     fn update(&mut self, p_in: &Vec<&[u8]>, channel_index: usize, add: bool) {
         let mut inc = 1;
         if !add {
@@ -96,6 +121,8 @@ impl PartialHistograms {
         }
     }
 
+    // Returns the number of pixels with a value of key in the kernel for the pixel at the given
+    // index
     fn get_count(&self, key: usize, index: usize) -> i32 {
         let mut count = self.data[self.n_half][key as usize];
         if index != self.n_half {
@@ -110,6 +137,20 @@ impl PartialHistograms {
 // Median filter functions
 ////////////////////////////
 
+/*
+ * The MedianHist struct:
+ *
+ * In addition to containing the partial histograms, this struct keeps track of each previous
+ * median of the kernel for each pixel. These medians are used as "pivots" to find the next
+ * median: to find the median of the kernel for a given pixel, instead of scanning its histogram
+ * starting from 0, we start from the median of the kernel of the previous pixel in that column.
+ * This value is typically much closer to the current median since the majority of the pixels
+ * in the previous and current kernels are the same, which makes scanning the histogram much
+ * quicker. The "sums", or the number of values in the current histogram that are less than the
+ * previous median, is used to determine if the current median is greater than or less than the
+ * previous median, thus determining if the current histogram should be scanned upwards or
+ * downwards from the previous median, respectively, to find the current median.
+ */
 #[derive(Debug, Clone)]
 struct MedianHist {
     data: PartialHistograms,
@@ -158,7 +199,7 @@ impl MedianHist {
             inc *= -1;
         }
 
-        // Update sums
+        // Update the number of values less than the previous median
         if !self.pivots.is_empty() {
             for n in 0..self.data.n_cols {
                 for i in n..(n + self.data.size) {
@@ -173,7 +214,9 @@ impl MedianHist {
 
 fn process_cols_med(input: &Image<u8>, output: &mut Image<u8>, radius: u32, n_cols: usize, x: u32) {
     let size = 2 * radius + 1;
-    let center = ((size * size) / 2 + 1) as i32;
+    let center = ((size * size) / 2 + 1) as i32; // Half the number of pixels in a kernel. If
+                                                      // all the pixels in the kernel were sorted,
+                                                      // the index of the median would be (center - 1).
     let (width, height, channels) = input.info().whc();
     let mut histograms = vec![MedianHist::new(radius as usize, n_cols); channels as usize];
 
@@ -252,10 +295,14 @@ fn process_row_med(output: &mut Image<u8>, histograms: &mut Vec<MedianHist>, cen
     for i in 0..n_cols {
         let mut p_out = Vec::with_capacity(channels);
         for c in 0..channels {
-            let pivot = histograms[c].pivots()[i];
-            let mut sum = histograms[c].sums()[i];
-
-            if sum < center {
+            let pivot = histograms[c].pivots()[i]; // Get the previous median
+            let mut sum = histograms[c].sums()[i]; // Get the number of values less than
+                                                        // the previous median
+
+            if sum == center { // The current median is equal to the previous median
+                p_out.push(pivot);
+            } else if sum < center { // The current median is greater than the previous median,
+                                     // so the histogram should be scanned upwards
                 for key in pivot..=255 {
                     let add = histograms[c].data().get_count(key as usize, i);
 
@@ -267,7 +314,8 @@ fn process_row_med(output: &mut Image<u8>, histograms: &mut Vec<MedianHist>, cen
 
                     sum += add;
                 }
-            } else {
+            } else { // The current median is less than the previous median, so the histogram
+                     // should be scanned downwards
                 for key in (0..pivot).rev() {
                     sum -= histograms[c].data().get_count(key as usize, i);
 
@@ -312,11 +360,14 @@ fn set_pivots_med(histograms: &mut Vec<MedianHist>, pivots: &Vec<u8>, index: usi
 #[derive(Debug, Clone)]
 struct MeanHist {
     data: PartialHistograms,
-    sums: Vec<i32>,
-    lower: Vec<Vec<u8>>,
-    upper: Vec<Vec<u8>>,
-    trim: usize,
-    len: f32,
+    sums: Vec<i32>, // The sum of all the participating pixel values in the kernel for each pixel
+    lower: Vec<Vec<u8>>, // Vectors of all the lowest discarded pixel values in the kernel
+                         // for each pixel
+    upper: Vec<Vec<u8>>, // Vectors of all the highest discarded pixel values in the kernel for
+                         // each pixel
+    trim: usize, // The number of pixel values discarded at the low and high ends of each kernel
+                 // (equal to half of alpha)
+    len: f32, // The number of participating pixel values in each kernel
 }
 
 impl MeanHist {
@@ -344,6 +395,7 @@ impl MeanHist {
         self.upper = vec![Vec::with_capacity(self.trim); self.data.n_cols];
     }
 
+    // By some miracle, this seems to work!
     fn update(&mut self, p_in: &Vec<&[u8]>, channel_index: usize, add: bool) {
         if !self.sums.is_empty() {
             if add {
diff --git a/src/filter/mod.rs b/src/filter/mod.rs
index 1785750..5971cf8 100644
--- a/src/filter/mod.rs
+++ b/src/filter/mod.rs
@@ -1,16 +1,18 @@
 //! A module for image filtering operations
 
-pub use self::median::*;
 pub use self::bilateral::*;
+pub use self::edge::*;
+pub use self::median::*;
 
 mod median;
 mod bilateral;
+mod edge;
 
-use crate::{util, colorspace, error, math};
-use crate::image::{Image, BaseImage, Number};
-use crate::util::constants::{K_SOBEL_1D_VERT, K_SOBEL_1D_HORZ, K_UNSHARP_MASKING, K_SHARPEN, K_PREWITT_1D_VERT, K_PREWITT_1D_HORZ};
+use crate::{error, util};
 use crate::enums::Thresh;
 use crate::error::ImgProcResult;
+use crate::image::{BaseImage, Image, Number};
+use crate::util::constants::{K_SHARPEN, K_UNSHARP_MASKING};
 
 #[cfg(feature = "rayon")]
 use rayon::prelude::*;
@@ -30,7 +32,7 @@ pub fn filter_1d(input: &Image<f64>, kernel: &[f64], is_vert: bool) -> ImgProcRe
 
     for y in 0..height {
         for x in 0..width {
-            let pixel = math::apply_1d_kernel(input.get_neighborhood_1d(x, y,
+            let pixel = util::apply_1d_kernel(input.get_neighborhood_1d(x, y,
                                                         kernel.len() as u32, is_vert), kernel)?;
             output.set_pixel(x, y, &pixel);
         }
@@ -51,7 +53,7 @@ pub fn filter_1d(input: &Image<f64>, kernel: &[f64], is_vert: bool) -> ImgProcRe
         .into_par_iter()
         .map(|i| {
             let (x, y) = util::get_2d_coords(i, width);
-            math::apply_1d_kernel(input.get_neighborhood_1d(x, y,kernel.len() as u32, is_vert), kernel).unwrap()
+            util::apply_1d_kernel(input.get_neighborhood_1d(x, y,kernel.len() as u32, is_vert), kernel).unwrap()
         })
         .collect();
 
@@ -80,7 +82,7 @@ pub fn unseparable_filter(input: &Image<f64>, kernel: &[f64]) -> ImgProcResult<I
 
     for y in 0..height {
         for x in 0..width {
-            let pixel = math::apply_2d_kernel(input.get_neighborhood_2d(x, y, size), kernel)?;
+            let pixel = util::apply_2d_kernel(input.get_neighborhood_2d(x, y, size), kernel)?;
             output.set_pixel(x, y, &pixel);
         }
     }
@@ -101,7 +103,7 @@ pub fn unseparable_filter(input: &Image<f64>, kernel: &[f64]) -> ImgProcResult<I
         .into_par_iter()
         .map(|i| {
             let (x, y) = util::get_2d_coords(i, width);
-            math::apply_2d_kernel(input.get_neighborhood_2d(x, y, size), kernel).unwrap()
+            util::apply_2d_kernel(input.get_neighborhood_2d(x, y, size), kernel).unwrap()
         })
         .collect();
 
@@ -113,7 +115,7 @@ pub fn linear_filter(input: &Image<f64>, kernel: &[f64]) -> ImgProcResult<Image<
     error::check_odd(kernel.len(), "kernel length")?;
     error::check_square(kernel.len() as f64, "kernel length")?;
 
-    let separable = math::separate_kernel(kernel);
+    let separable = util::separate_kernel(kernel);
     match separable {
         Some((vert, horz)) => Ok(separable_filter(input, &vert, &horz)?),
         None => Ok(unseparable_filter(input, &kernel)?)
@@ -147,8 +149,8 @@ pub fn weighted_avg_filter(input: &Image<f64>, size: u32, weight: u32) -> ImgPro
 }
 
 /// Applies a Gaussian blur using a `size x size` kernel
-pub fn gaussian_blur(input: &Image<f64>, size: u32, std_dev: f64) -> ImgProcResult<Image<f64>> {
-    let kernel = util::generate_gaussian_kernel(size, std_dev)?;
+pub fn gaussian_blur(input: &Image<f64>, size: u32, sigma: f64) -> ImgProcResult<Image<f64>> {
+    let kernel = util::generate_gaussian_kernel(size, sigma)?;
     Ok(linear_filter(input, &kernel)?)
 }
 
@@ -166,48 +168,13 @@ pub fn unsharp_masking(input: &Image<f64>) -> ImgProcResult<Image<f64>> {
     Ok(unseparable_filter(input, &K_UNSHARP_MASKING)?)
 }
 
-////////////////////
-// Edge detection
-////////////////////
-
-/// Applies a separable derivative mask; first converts `input` to grayscale
-pub fn derivative_mask(input: &Image<f64>, vert_kernel: &[f64], horz_kernel: &[f64]) -> ImgProcResult<Image<f64>> {
-    let gray = colorspace::rgb_to_grayscale_f64(input);
-    let img_x = separable_filter(&gray, &vert_kernel, &horz_kernel)?;
-    let img_y = separable_filter(&gray, &horz_kernel, &vert_kernel)?;
-
-    let mut output = Image::blank(gray.info());
-
-    for i in 0..(output.info().full_size() as usize) {
-        output.set_pixel_indexed(i, &[(img_x[i][0].powf(2.0) + img_y[i][0].powf(2.0)).sqrt()]);
-    }
-
-    Ok(output)
-}
-
-/// Applies the Prewitt operator
-pub fn prewitt(input: &Image<f64>) -> ImgProcResult<Image<f64>> {
-    Ok(derivative_mask(input, &K_PREWITT_1D_VERT, &K_PREWITT_1D_HORZ)?)
-}
-
-/// Applies the Sobel operator
-pub fn sobel(input: &Image<f64>) -> ImgProcResult<Image<f64>> {
-    Ok(derivative_mask(input, &K_SOBEL_1D_VERT, &K_SOBEL_1D_HORZ)?)
-}
-
-/// Applies a Sobel operator with weight `weight`
-pub fn sobel_weighted(input: &Image<f64>, weight: u32) -> ImgProcResult<Image<f64>> {
-    let vert_kernel = vec![1.0, weight as f64, 1.0];
-    Ok(derivative_mask(input, &vert_kernel, &K_SOBEL_1D_HORZ)?)
-}
-
 //////////////////
 // Thresholding
 //////////////////
 
 /// Performs a thresholding operation based on `method`
 pub fn threshold(input: &Image<f64>, threshold: f64, max: f64, method: Thresh) -> ImgProcResult<Image<f64>> {
-    error::check_grayscale(input.info().channels, input.info().alpha)?;
+    error::check_grayscale(input)?;
 
     match method {
         Thresh::Binary => {
diff --git a/src/image/pixel_iter.rs b/src/image/pixel_iter.rs
index af83648..6dba5e5 100644
--- a/src/image/pixel_iter.rs
+++ b/src/image/pixel_iter.rs
@@ -1,6 +1,30 @@
 use crate::image::{Image, BaseImage, Number};
 
-/// A struct representing a pixel iterator for an image
+/// A struct representing a pixel iterator for an image. `next()` returns a tuple containing the
+/// x-coordinate, y-coordinate, and a slice representing the pixel at that coordinate, in that
+/// order.
+///
+/// # Examples
+/// ```rust
+/// # fn main() {
+/// use imgproc_rs::image::{Image, BaseImage};
+///
+/// // Create an image
+/// let img = Image::from_vec(2, 2, 3, false, vec![1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]);
+///
+/// // Print pixels with corresponding coordinates using the pixel iterator
+/// for vals in img.into_iter() {
+///     print!("(x: {}, y: {}), pixel: (", vals.0, vals.1);
+///
+///     for i in 0..(img.info().channels as usize) {
+///         print!("{}, ", vals.2[i]);
+///     }
+///
+///     print!(")");
+///     println!();
+/// }
+/// # }
+/// ```
 #[derive(Debug, Clone)]
 pub struct PixelIter<'a, T: Number> {
     image: &'a Image<T>,
diff --git a/src/lib.rs b/src/lib.rs
index c19cac8..5c89ff7 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -4,7 +4,6 @@ pub mod io;
 pub mod error;
 pub mod util;
 pub mod enums;
-pub mod math;
 pub mod colorspace;
 pub mod tone;
 pub mod filter;
diff --git a/src/morphology.rs b/src/morphology.rs
index 800fc89..5d5eb52 100644
--- a/src/morphology.rs
+++ b/src/morphology.rs
@@ -5,12 +5,12 @@ use crate::image::{Image, BaseImage};
 /// Erodes a binary image (grayscale image with pixel values of 0 or 255) using a kernel of size
 /// `(2 * radius + 1) x (2 * radius + 1)`
 pub fn erode(input: &Image<u8>, radius: u32) -> ImgProcResult<Image<u8>> {
-    error::check_grayscale(input.info().channels, input.info().alpha)?;
+    error::check_grayscale(input)?;
 
     let (width, height) =  input.info().wh();
     let size = 2 * radius + 1;
     let max_sum = (size * size * 255) as f64;
-    let table = util::summed_area_table(&input.clone().into());
+    let table = util::generate_summed_area_table(&input.clone().into());
     let mut output = Image::blank(input.info());
 
     for y in 0..height {
@@ -45,10 +45,10 @@ pub fn erode(input: &Image<u8>, radius: u32) -> ImgProcResult<Image<u8>> {
 /// Dilates a binary image (grayscale image with pixel values of 0 or 255) using a kernel of size
 /// `(2 * radius + 1) x (2 * radius + 1)`
 pub fn dilate(input: &Image<u8>, radius: u32) -> ImgProcResult<Image<u8>> {
-    error::check_grayscale(input.info().channels, input.info().alpha)?;
+    error::check_grayscale(input)?;
 
     let (width, height) =  input.info().wh();
-    let table = util::summed_area_table(&input.clone().into());
+    let table = util::generate_summed_area_table(&input.clone().into());
     let mut output = Image::blank(input.info());
 
     for y in 0..height {
@@ -83,10 +83,10 @@ pub fn dilate(input: &Image<u8>, radius: u32) -> ImgProcResult<Image<u8>> {
 /// Sets output pixel to the majority-valued pixel in the input image under a kernel of size
 /// `(2 * radius + 1) x (2 * radius + 1)`
 pub fn majority(input: &Image<u8>, radius: u32) -> ImgProcResult<Image<u8>> {
-    error::check_grayscale(input.info().channels, input.info().alpha)?;
+    error::check_grayscale(input)?;
 
     let (width, height) =  input.info().wh();
-    let table = util::summed_area_table(&input.clone().into());
+    let table = util::generate_summed_area_table(&input.clone().into());
     let mut output = Image::blank(input.info());
 
     for y in 0..height {
@@ -131,12 +131,12 @@ pub fn close(input: &Image<u8>, radius: u32) -> ImgProcResult<Image<u8>> {
 /// Returns the difference between dilation and erosion of the image
 #[allow(unused_parens)]
 pub fn gradient(input: &Image<u8>, radius: u32) -> ImgProcResult<Image<u8>> {
-    error::check_grayscale(input.info().channels, input.info().alpha)?;
+    error::check_grayscale(input)?;
 
     let (width, height) =  input.info().wh();
     let size = 2 * radius + 1;
     let max_sum = (size * size * 255) as f64;
-    let table = util::summed_area_table(&input.clone().into());
+    let table = util::generate_summed_area_table(&input.clone().into());
     let mut output = Image::blank(input.info());
 
     for y in 0..height {
diff --git a/src/tone.rs b/src/tone.rs
index 799fad3..1899444 100644
--- a/src/tone.rs
+++ b/src/tone.rs
@@ -3,21 +3,23 @@
 use crate::{util, colorspace, error};
 use crate::enums::{Tone, White};
 use crate::image::Image;
-use crate::error::{ImgProcError, ImgProcResult};
+use crate::error::ImgProcResult;
 
 use std::collections::HashMap;
 
-/// Adjusts brightness by adding `bias` to each RGB channel if `method` is `Tone::Rgb` or adding
-/// `bias` to the luminance value (Y) of `input` in CIE XYZ if `method` is `Tone::Xyz`
+/// Adjusts brightness by adding `bias` to each RGB channel if `method` is `Tone::Rgb`, or adding
+/// `bias` to the L* channel of `input` in CIELAB if `method` is `Tone::Lab`
+///
+/// # Arguments
+///
+/// * `bias` - Must be between 0 and 255 (inclusive)
 pub fn brightness(input: &Image<u8>, bias: i32, method: Tone) -> ImgProcResult<Image<u8>> {
-    if bias < 0 || bias > 255 {
-        return Err(ImgProcError::InvalidArgError("bias is not in range 0 to 255".to_string()));
-    }
+    error::check_in_range(bias, 0, 255, "bias")?;
 
     match method {
         Tone::Rgb => {
             let mut lookup_table: [u8; 256] = [0; 256];
-            util::create_lookup_table(&mut lookup_table, |i| {
+            util::generate_lookup_table(&mut lookup_table, |i| {
                 (i as i32 + bias).clamp(0, 255) as u8
             });
 
@@ -31,16 +33,19 @@ pub fn brightness(input: &Image<u8>, bias: i32, method: Tone) -> ImgProcResult<I
     }
 }
 
-/// Adjusts contrast by multiplying each RGB channel by `gain` if `method` is `Tone::Rgb` or
-/// multiplying the luminance value (Y) of `input` in CIE XYZ by `gain` if `method` is `Tone::Xyz`
-// gain > 0
+/// Adjusts contrast by multiplying each RGB channel by `gain` if `method` is `Tone::Rgb`, or
+/// multiplying the L* channel of `input` in CIELAB by `gain` if `method` is `Tone::Lab`
+///
+/// # Arguments
+///
+/// * `gain` - Must be between 0 and 1 (inclusive)
 pub fn contrast(input: &Image<u8>, gain: f64, method: Tone) -> ImgProcResult<Image<u8>> {
     error::check_non_neg(gain, "gain")?;
 
     match method {
         Tone::Rgb => {
             let mut lookup_table: [u8; 256] = [0; 256];
-            util::create_lookup_table(&mut lookup_table, |i| {
+            util::generate_lookup_table(&mut lookup_table, |i| {
                 (i as f64 * gain).round().clamp(0.0, 255.0) as u8
             });
 
@@ -55,7 +60,13 @@ pub fn contrast(input: &Image<u8>, gain: f64, method: Tone) -> ImgProcResult<Ima
 }
 
 /// Adjusts saturation by adding `saturation` to the saturation value (S) of `input` in HSV
+///
+/// # Arguments
+///
+/// * `saturation` - Must be between 0 and 255 (inclusive)
 pub fn saturation(input: &Image<u8>, saturation: i32) -> ImgProcResult<Image<u8>> {
+    error::check_in_range(saturation, 0, 255, "saturation")?;
+
     let mut hsv = colorspace::rgb_to_hsv(input);
     hsv.edit_channel(|s| (s + (saturation as f64 / 255.0)) as f64, 1);
 
@@ -63,6 +74,10 @@ pub fn saturation(input: &Image<u8>, saturation: i32) -> ImgProcResult<Image<u8>
 }
 
 /// Performs a gamma correction. `max` indicates the maximum allowed pixel value of the image
+///
+/// # Arguments
+///
+/// * `gamma` - Must be non-negative
 pub fn gamma(input: &Image<u8>, gamma: f64, max: u8) -> ImgProcResult<Image<u8>> {
     error::check_non_neg(gamma, "gamma")?;
 
@@ -78,12 +93,12 @@ pub fn gamma(input: &Image<u8>, gamma: f64, max: u8) -> ImgProcResult<Image<u8>>
 /// * `alpha` - Represents the amount of equalization, where 0 corresponds to no equalization and
 /// 1 corresponds to full equalization
 /// * `ref_white` - An enum representing the reference white value of the image
-/// * `precision` - See the function `util::generate_histogram_percentiles`
+/// * `precision` - Must be non-negative. See
+/// [`generate_histogram_percentiles`](../util/fn.generate_histogram_percentiles.html) for a
+/// complete description
 pub fn histogram_equalization(input: &Image<u8>, alpha: f64, ref_white: &White, precision: f64) -> ImgProcResult<Image<u8>> {
     error::check_non_neg(precision, "precision")?;
-    if alpha < 0.0 || alpha > 1.0 {
-        return Err(ImgProcError::InvalidArgError("alpha is not in range 0 to 1".to_string()));
-    }
+    error::check_in_range(alpha, 0.0, 1.0, "alpha")?;
 
     let mut lab = colorspace::srgb_to_lab(input, ref_white);
     let mut percentiles = HashMap::new();
diff --git a/src/transform.rs b/src/transform.rs
index 6327f3d..17c2881 100644
--- a/src/transform.rs
+++ b/src/transform.rs
@@ -1,15 +1,14 @@
 //! A module for image transformation operations
 
-use crate::{math, error};
-#[cfg(feature = "rayon")]
-use crate::util;
-use crate::image::{Number, Image, ImageInfo, BaseImage};
-use crate::error::{ImgProcResult, ImgProcError};
-use crate::enums::{Scale, Refl};
-
 #[cfg(feature = "rayon")]
 use rayon::prelude::*;
 
+use crate::enums::{Refl, Scale};
+use crate::error;
+use crate::error::{ImgProcError, ImgProcResult};
+use crate::image::{BaseImage, Image, ImageInfo, Number};
+use crate::util;
+
 /// Crops an image to a rectangle with upper left corner located at `(x, y)` with width `width`
 /// and height `height`
 #[cfg(not(feature = "rayon"))]
@@ -221,15 +220,15 @@ pub fn rotate(input: &Image<f64>, degrees: f64) -> ImgProcResult<Image<f64>> {
     let mat = [cos, -sin, sin, cos];
 
     // Compute dimensions of output image
-    let coords1 = math::vector_mul(&mat, &[-(x as f64), y as f64])?;
-    let coords2 = math::vector_mul(&mat, &[(w_in - x) as f64, y as f64])?;
-    let coords3 = math::vector_mul(&mat, &[-(x as f64), (y as f64) - (h_in as f64)])?;
-    let coords4 = math::vector_mul(&mat, &[(w_in - x) as f64, (y as f64) - (h_in as f64)])?;
+    let coords1 = util::vector_mul(&mat, &[-(x as f64), y as f64])?;
+    let coords2 = util::vector_mul(&mat, &[(w_in - x) as f64, y as f64])?;
+    let coords3 = util::vector_mul(&mat, &[-(x as f64), (y as f64) - (h_in as f64)])?;
+    let coords4 = util::vector_mul(&mat, &[(w_in - x) as f64, (y as f64) - (h_in as f64)])?;
 
-    let x_max = math::max_4(coords1[0], coords2[0], coords3[0], coords4[0]);
-    let x_min = math::min_4(coords1[0], coords2[0], coords3[0], coords4[0]);
-    let y_max = math::max_4(coords1[1], coords2[1], coords3[1], coords4[1]);
-    let y_min = math::min_4(coords1[1], coords2[1], coords3[1], coords4[1]);
+    let x_max = util::max_4(coords1[0], coords2[0], coords3[0], coords4[0]);
+    let x_min = util::min_4(coords1[0], coords2[0], coords3[0], coords4[0]);
+    let y_max = util::max_4(coords1[1], coords2[1], coords3[1], coords4[1]);
+    let y_min = util::min_4(coords1[1], coords2[1], coords3[1], coords4[1]);
 
     let w_out = (x_max - x_min) as u32;
     let h_out = (y_max - y_min) as u32;
@@ -242,7 +241,7 @@ pub fn rotate(input: &Image<f64>, degrees: f64) -> ImgProcResult<Image<f64>> {
             let x1 = (i as f64) - (x as f64);
             let y1 = (y as f64) - (j as f64);
 
-            let mut coords = math::vector_mul(&mat, &[x1, y1])?;
+            let mut coords = util::vector_mul(&mat, &[x1, y1])?;
 
             coords[0] += x_max - 1.0;
             coords[1] = y_max - coords[1] - 1.0;
@@ -302,7 +301,7 @@ pub fn shear(input: &Image<f64>, shear_x: f64, shear_y: f64) -> ImgProcResult<Im
 
     for y in 0..h_in {
         for x in 0..w_in {
-            let mut coords = math::vector_mul(&mat, &[x as f64, y as f64])?;
+            let mut coords = util::vector_mul(&mat, &[x as f64, y as f64])?;
 
             if shear_x > 0.0 {
                 coords[0] += offset_x;
@@ -473,8 +472,8 @@ fn interpolate_bicubic(input: &Image<f64>, x_factor: f64, y_factor: f64, x: u32,
         for n in -1..3 {
             let y_clamp = (y_in + (n as f64)).clamp(0.0, input.info().height as f64 - 1.0) as u32;
             let p_in = input.get_pixel_unchecked(x_clamp, y_clamp);
-            let r = math::cubic_weighting_fn((m as f64) - delta_x)
-                * math::cubic_weighting_fn(delta_y - (n as f64));
+            let r = util::cubic_weighting_fn((m as f64) - delta_x)
+                * util::cubic_weighting_fn(delta_y - (n as f64));
 
             for c in 0..(input.info().channels as usize) {
                 p_out[c] += p_in[c] * r;
@@ -498,8 +497,8 @@ fn interpolate_lanczos(input: &Image<f64>, x_factor: f64, y_factor: f64, size: u
         for j in (1 - (size as i32))..(size as i32 + 1) {
             let y_clamp = (y_in + (j as f64)).clamp(0.0, input.info().height as f64 - 1.0) as u32;
             let p_in = input.get_pixel_unchecked(x_clamp, y_clamp);
-            let lanczos = math::lanczos_kernel(delta_x - (i as f64), size as f64)
-                * math::lanczos_kernel(delta_y - (j as f64), size as f64);
+            let lanczos = util::lanczos_kernel(delta_x - (i as f64), size as f64)
+                * util::lanczos_kernel(delta_y - (j as f64), size as f64);
 
             for c in 0..(input.info().channels as usize) {
                 p_out[c] += p_in[c] * lanczos;
diff --git a/src/util/constants.rs b/src/util/constants.rs
index 2cf71d2..0ae61f8 100644
--- a/src/util/constants.rs
+++ b/src/util/constants.rs
@@ -62,4 +62,7 @@ pub const K_SOBEL_1D_HORZ: [f64; 3] = [-1.0, 0.0, 1.0];
 pub const K_PREWITT_1D_VERT: [f64; 3] = [1.0, 1.0, 1.0];
 
 /// 1D horizontal kernel for the Prewitt operator
-pub const K_PREWITT_1D_HORZ: [f64; 3] = [-1.0, 0.0, 1.0];
\ No newline at end of file
+pub const K_PREWITT_1D_HORZ: [f64; 3] = [-1.0, 0.0, 1.0];
+
+/// Laplacian discrete approximation kernel
+pub const K_LAPLACIAN: [f64; 9] = [0.0, -1.0, 0.0, -1.0, 4.0, -1.0, 0.0, -1.0, 0.0];
\ No newline at end of file
diff --git a/src/math.rs b/src/util/math.rs
similarity index 90%
rename from src/math.rs
rename to src/util/math.rs
index 3736c95..064e7ae 100644
--- a/src/math.rs
+++ b/src/util/math.rs
@@ -1,5 +1,3 @@
-//! A module for utility math functions
-
 use crate::error;
 use crate::image::{Number, SubImage, BaseImage};
 use crate::error::ImgProcResult;
@@ -177,6 +175,28 @@ pub fn cubic_weighting_fn(x: f64) -> f64 {
         - 4.0 * clamp_zero(x - 1.0).powf(3.0))
 }
 
+/// A helper function for the colorspace conversion from CIE XYZ to CIELAB
+pub fn xyz_to_lab_fn(num: f64) -> f64 {
+    let d: f64 = 6.0 / 29.0;
+
+    if num > d.powf(3.0) {
+        num.powf(1.0 / 3.0)
+    } else {
+        (num / (3.0 * d * d)) + (4.0 / 29.0)
+    }
+}
+
+/// A helper function for the colorspace conversion from CIELAB to CIE XYZ
+pub fn lab_to_xyz_fn(num: f64) -> f64 {
+    let d: f64 = 6.0 / 29.0;
+
+    if num > d {
+        num.powf(3.0)
+    } else {
+        3.0 * d * d * (num - (4.0 / 29.0))
+    }
+}
+
 /// Returns 0 if `x` is less than 0; `x` if not
 pub fn clamp_zero(x: f64) -> f64 {
     if x <= 0.0 {
diff --git a/src/util/mod.rs b/src/util/mod.rs
index 75dd831..70e09c1 100644
--- a/src/util/mod.rs
+++ b/src/util/mod.rs
@@ -1,49 +1,31 @@
 //! A module for image utility functions
 
-pub mod constants;
+pub use self::math::*;
+
+mod math;
+
+use std::collections::{BTreeMap, HashMap};
+use std::f64::consts::{E, PI};
 
-use crate::{error, math};
 use crate::enums::White;
+use crate::error;
 use crate::error::ImgProcResult;
 use crate::image::{BaseImage, Image, Number};
 
-use std::collections::{BTreeMap, HashMap};
-use std::f64::consts::{E, PI};
+pub mod constants;
 
 ////////////////////////////
 // Image helper functions
 ////////////////////////////
 
 /// Returns a tuple representing the XYZ tristimulus values for a given reference white value
-pub fn generate_xyz_tristimulus_vals(ref_white: &White) -> (f64, f64, f64) {
+pub fn xyz_tristimulus_vals(ref_white: &White) -> (f64, f64, f64) {
     return match ref_white {
         White::D50 => (96.4212, 100.0, 82.5188),
         White::D65 => (95.0489, 100.0, 108.8840),
     }
 }
 
-/// A helper function for the colorspace conversion from CIE XYZ to CIELAB
-pub fn xyz_to_lab_fn(num: f64) -> f64 {
-    let d: f64 = 6.0 / 29.0;
-
-    if num > d.powf(3.0) {
-        num.powf(1.0 / 3.0)
-    } else {
-        (num / (3.0 * d * d)) + (4.0 / 29.0)
-    }
-}
-
-/// A helper function for the colorspace conversion from CIELAB to CIE XYZ
-pub fn lab_to_xyz_fn(num: f64) -> f64 {
-    let d: f64 = 6.0 / 29.0;
-
-    if num > d {
-        num.powf(3.0)
-    } else {
-        3.0 * d * d * (num - (4.0 / 29.0))
-    }
-}
-
 /// A helper function for histogram equalization
 ///
 /// # Arguments
@@ -73,7 +55,7 @@ pub fn generate_histogram_percentiles(input: &Image<f64>, percentiles: &mut Hash
 }
 
 /// Populates `table` with the appropriate values based on function `f`
-pub fn create_lookup_table<T: Number, F>(table: &mut [T; 256], f: F)
+pub fn generate_lookup_table<T: Number, F>(table: &mut [T; 256], f: F)
     where F: Fn(u8) -> T {
     for i in 0..256 {
         table[i] = f(i as u8);
@@ -81,7 +63,31 @@ pub fn create_lookup_table<T: Number, F>(table: &mut [T; 256], f: F)
 }
 
 /// Generates a Gaussian kernel
-pub fn generate_gaussian_kernel(size: u32, std_dev: f64) -> ImgProcResult<Vec<f64>> {
+pub fn generate_gaussian_kernel(size: u32, sigma: f64) -> ImgProcResult<Vec<f64>> {
+    error::check_odd(size, "size")?;
+
+    let mut filter = vec![0.0; (size * size) as usize];
+    let k = ((size - 1) / 2) as i32;
+
+    for i in 0..(size as i32) {
+        for j in 0..(size as i32) {
+            if i <= j {
+                let num = (1.0 / (2.0 * PI * sigma * sigma)) *
+                    (E.powf(-(((i - k) * (i - k) + (j - k) * (j - k)) as f64) / (2.0 * sigma * sigma)));
+                filter[(i * size as i32 + j) as usize] = num;
+
+                if i != j {
+                    filter[(j * size as i32 + i) as usize] = num;
+                }
+            }
+        }
+    }
+
+    Ok(filter)
+}
+
+/// Generates a Laplacian of Gaussian kernel
+pub fn generate_log_kernel(size: u32, sigma: f64) -> ImgProcResult<Vec<f64>> {
     error::check_odd(size, "size")?;
 
     let mut filter = vec![0.0; (size * size) as usize];
@@ -90,8 +96,8 @@ pub fn generate_gaussian_kernel(size: u32, std_dev: f64) -> ImgProcResult<Vec<f6
     for i in 0..(size as i32) {
         for j in 0..(size as i32) {
             if i <= j {
-                let num = (1.0 / (2.0 * PI * std_dev * std_dev)) *
-                    (E.powf(-(((i - k) * (i - k) + (j - k) * (j - k)) as f64) / (2.0 * std_dev * std_dev)));
+                let exp = -(((i - k) * (i - k) + (j - k) * (j - k)) as f64) / (2.0 * sigma * sigma);
+                let num = (-1.0 / (PI * sigma.powf(4.0))) * (1.0 - exp) * (E.powf(exp));
                 filter[(i * size as i32 + j) as usize] = num;
 
                 if i != j {
@@ -154,7 +160,7 @@ pub fn generate_spatial_mat(size: u32, spatial: f64) -> ImgProcResult<Vec<f64>>
 
 /// Generates a summed-area table in the format of another `Image` of the same type and dimensions
 /// as `input`
-pub fn summed_area_table(input: &Image<f64>) -> Image<f64> {
+pub fn generate_summed_area_table(input: &Image<f64>) -> Image<f64> {
     let mut output = Image::blank(input.info());
     let (width, height, channels) = input.info().whc();
     let zeroes = vec![0.0; channels as usize];
diff --git a/tests/colorspace_test.rs b/tests/colorspace_test.rs
index 55b4145..7b87fe9 100644
--- a/tests/colorspace_test.rs
+++ b/tests/colorspace_test.rs
@@ -42,7 +42,7 @@ fn srgb_to_xyz_test() {
     let proc = colorspace::srgb_to_xyz(&img);
     println!("processing: {}", now.elapsed().unwrap().as_millis());
 
-    write(&convert::scale_channels(&proc, 1.0, 255.0).unwrap().into(), "images/tests/colorspace/srgb_xyz.png").unwrap();
+    write(&convert::scale_channels(&proc, 0.0, 0.0, 1.0, 255.0).unwrap().into(), "images/tests/colorspace/srgb_xyz.png").unwrap();
 }
 
 // #[test]
@@ -50,7 +50,7 @@ fn xyz_to_srgb_test() {
     let img: Image<f64> = setup("images/tests/colorspace/srgb_xyz.png").unwrap().into();
 
     let now = SystemTime::now();
-    let proc = colorspace::xyz_to_srgb(&convert::scale_channels(&img, 255.0, 1.0).unwrap());
+    let proc = colorspace::xyz_to_srgb(&convert::scale_channels(&img, 0.0, 0.0, 255.0, 1.0).unwrap());
     println!("processing: {}", now.elapsed().unwrap().as_millis());
 
     write(&proc, "images/tests/colorspace/xyz_srgb.png").unwrap();
@@ -94,7 +94,7 @@ fn rgb_to_hsv_test() {
     let proc = colorspace::rgb_to_hsv(&img);
     println!("processing: {}", now.elapsed().unwrap().as_millis());
 
-    write(&convert::scale_channels(&proc, 1.0, 255.0).unwrap().into(), "images/tests/colorspace/rgb_hsv.png").unwrap();
+    write(&convert::scale_channels(&proc, 0.0, 0.0, 1.0, 255.0).unwrap().into(), "images/tests/colorspace/rgb_hsv.png").unwrap();
 }
 
 // #[test]
@@ -102,7 +102,7 @@ fn hsv_to_rgb_test() {
     let img: Image<f64> = setup("images/tests/colorspace/rgb_hsv.png").unwrap().into();
 
     let now = SystemTime::now();
-    let proc = colorspace::hsv_to_rgb(&convert::scale_channels(&img, 255.0, 1.0).unwrap());
+    let proc = colorspace::hsv_to_rgb(&convert::scale_channels(&img, 0.0, 0.0, 255.0, 1.0).unwrap());
     println!("processing: {}", now.elapsed().unwrap().as_millis());
 
     write(&proc, "images/tests/colorspace/hsv_rgb.png").unwrap();
diff --git a/tests/filter_test.rs b/tests/filter_test.rs
index d48499e..ad3a817 100644
--- a/tests/filter_test.rs
+++ b/tests/filter_test.rs
@@ -133,6 +133,28 @@ fn sobel_weighted() {
     write(&filtered.into(), "images/tests/filter/sobel_weighted.png").unwrap();
 }
 
+// #[test]
+fn laplacian() {
+    let img: Image<f64> = colorspace::rgb_to_grayscale(&setup("images/poppy.jpg").unwrap()).into();
+
+    let now = SystemTime::now();
+    let filtered = filter::laplacian(&img).unwrap();
+    println!("processing: {}", now.elapsed().unwrap().as_millis());
+
+    write(&filter::normalize_laplacian(&filtered).unwrap(), "images/tests/filter/laplacian.png").unwrap();
+}
+
+// #[test]
+fn laplacian_of_gaussian() {
+    let img: Image<f64> = colorspace::rgb_to_grayscale(&setup("images/scaled.png").unwrap()).into();
+
+    let now = SystemTime::now();
+    let filtered = filter::laplacian_of_gaussian(&img, 7, 1.0).unwrap();
+    println!("processing: {}", now.elapsed().unwrap().as_millis());
+
+    write(&filter::normalize_laplacian(&filtered).unwrap(), "images/tests/filter/laplacian_of_gaussian.png").unwrap();
+}
+
 // #[test]
 fn threshold_test() {
     let img: Image<f64> = colorspace::rgb_to_grayscale(&setup(PATH).unwrap()).into();
diff --git a/tests/math_test.rs b/tests/math_test.rs
index 9bbac57..55263a3 100644
--- a/tests/math_test.rs
+++ b/tests/math_test.rs
@@ -1,4 +1,4 @@
-use imgproc_rs::math;
+use imgproc_rs::util;
 use imgproc_rs::util::constants::K_GAUSSIAN_BLUR_2D_3;
 use imgproc_rs::image::SubImage;
 
@@ -6,7 +6,7 @@ use imgproc_rs::image::SubImage;
 fn vector_mul_test() {
     let mat = vec![1, 2, 3, 4, 5, 6, 7, 8, 9];
     let vec = vec![1, 2, 3];
-    let res = math::vector_mul(&mat, &vec).unwrap();
+    let res = util::vector_mul(&mat, &vec).unwrap();
 
     assert_eq!(vec![14, 32, 50], res);
 }
@@ -14,31 +14,31 @@ fn vector_mul_test() {
 #[test]
 fn max_test() {
     // Test max_3()
-    assert_eq!(3.0, math::max_3(1.0, 2.0, 3.0));
-    assert_eq!(1.0, math::max_3(1.0, 1.0, 1.0));
-    assert_eq!(3.0, math::max_3(1.0, 3.0, 3.0));
-    assert_eq!(3.0, math::max_3(1.0, 1.0, 3.0));
+    assert_eq!(3.0, util::max_3(1.0, 2.0, 3.0));
+    assert_eq!(1.0, util::max_3(1.0, 1.0, 1.0));
+    assert_eq!(3.0, util::max_3(1.0, 3.0, 3.0));
+    assert_eq!(3.0, util::max_3(1.0, 1.0, 3.0));
 
     // Test max_4()
-    assert_eq!(4.0, math::max_4(1.0, 2.0, 3.0, 4.0));
-    assert_eq!(1.0, math::max_4(1.0, 1.0, 1.0, 1.0));
-    assert_eq!(3.0, math::max_4(1.0, 2.0, 2.0, 3.0));
-    assert_eq!(3.0, math::max_4(1.0, 2.0, 3.0, 3.0));
+    assert_eq!(4.0, util::max_4(1.0, 2.0, 3.0, 4.0));
+    assert_eq!(1.0, util::max_4(1.0, 1.0, 1.0, 1.0));
+    assert_eq!(3.0, util::max_4(1.0, 2.0, 2.0, 3.0));
+    assert_eq!(3.0, util::max_4(1.0, 2.0, 3.0, 3.0));
 }
 
 #[test]
 fn min_test() {
     // Test min_3()
-    assert_eq!(1.0, math::min_3(1.0, 2.0, 3.0));
-    assert_eq!(1.0, math::min_3(1.0, 1.0, 1.0));
-    assert_eq!(1.0, math::min_3(1.0, 3.0, 3.0));
-    assert_eq!(1.0, math::min_3(1.0, 1.0, 3.0));
+    assert_eq!(1.0, util::min_3(1.0, 2.0, 3.0));
+    assert_eq!(1.0, util::min_3(1.0, 1.0, 1.0));
+    assert_eq!(1.0, util::min_3(1.0, 3.0, 3.0));
+    assert_eq!(1.0, util::min_3(1.0, 1.0, 3.0));
 
     // Test min_4()
-    assert_eq!(1.0, math::min_4(1.0, 2.0, 3.0, 4.0));
-    assert_eq!(1.0, math::min_4(1.0, 1.0, 1.0, 1.0));
-    assert_eq!(1.0, math::min_4(1.0, 2.0, 2.0, 3.0));
-    assert_eq!(1.0, math::min_4(1.0, 2.0, 3.0, 3.0));
+    assert_eq!(1.0, util::min_4(1.0, 2.0, 3.0, 4.0));
+    assert_eq!(1.0, util::min_4(1.0, 1.0, 1.0, 1.0));
+    assert_eq!(1.0, util::min_4(1.0, 2.0, 2.0, 3.0));
+    assert_eq!(1.0, util::min_4(1.0, 2.0, 3.0, 3.0));
 }
 
 #[test]
@@ -48,7 +48,7 @@ fn apply_1d_kernel_test() {
                                &[2.0, 3.0, 4.0]];
     let subimg = SubImage::new(3, 1, 3, false, pixels);
     let kernel = [1.0, 2.0, 1.0];
-    let res = math::apply_1d_kernel(subimg, &kernel).unwrap();
+    let res = util::apply_1d_kernel(subimg, &kernel).unwrap();
 
     assert_eq!(vec![11.0, 15.0, 19.0], res);
 }
@@ -65,7 +65,7 @@ fn apply_2d_kernel_test() {
                       &[3.0, 5.0, 7.0],
                       &[1.0, 3.0, 5.0]];
     let subimg = SubImage::new(3, 3, 3, false, pixels);
-    let res = math::apply_2d_kernel(subimg, &K_GAUSSIAN_BLUR_2D_3).unwrap();
+    let res = util::apply_2d_kernel(subimg, &K_GAUSSIAN_BLUR_2D_3).unwrap();
 
     assert_eq!(vec![3.5625, 3.8125, 4.0625], res);
 }
\ No newline at end of file
diff --git a/tests/util_test.rs b/tests/util_test.rs
index 704d018..327964c 100644
--- a/tests/util_test.rs
+++ b/tests/util_test.rs
@@ -7,7 +7,7 @@ fn summed_area_table_test() {
                                   &[31.0, 2.0, 4.0, 33.0, 5.0, 36.0,
                                12.0, 26.0, 9.0, 10.0, 29.0, 25.0,
                                13.0, 17.0, 21.0, 22.0, 20.0, 18.0]);
-    let output = util::summed_area_table(&input);
+    let output = util::generate_summed_area_table(&input);
     let output_table = [31.0, 33.0, 37.0, 70.0, 75.0, 111.0,
                                     43.0, 71.0, 84.0, 127.0, 161.0, 222.0,
                                     56.0, 101.0, 135.0, 200.0, 254.0, 333.0];