diff --git a/turbidity/main.go b/turbidity/main.go
index f5cc1874..ad9c2e3d 100644
--- a/turbidity/main.go
+++ b/turbidity/main.go
@@ -12,58 +12,69 @@ import (
 )
 
 func main() {
-	n := 10
+	n_images := 6
+	n_samples := 10
 
 	// Load template
 	template := gocv.IMRead("template.jpg", gocv.IMReadGrayScale)
+	standard := gocv.IMRead("standard.jpg", gocv.IMReadGrayScale)
 
-	// Read images
-	images := make([]gocv.Mat, n)
+	images := make([][]gocv.Mat, n_images)
 
-	for i := 0; i < n; i++ {
-		images[i] = gocv.IMRead(fmt.Sprintf("images/generation-5/000%v.jpg", i), gocv.IMReadGrayScale)
+	//Load images
+	for i := 0; i < n_images; i++ {
+		images[i] = make([]gocv.Mat, n_samples)
+		for j := 0; j < n_samples; j++ {
+			images[i][j] = gocv.IMRead(fmt.Sprintf("images/t-%v/000%v.jpg", i, j), gocv.IMReadGrayScale)
+		}
 	}
 
-	// Calulate turbidity
-	ts := TurbiditySensor{template: template, k1: 90, k2: 90, scale: 5.0, alpha: 1.0}
+	// Create turbidity sensor
+	ts := TurbiditySensor{template: template, standard: standard, k1: 90, k2: 90, scale: 5.0, alpha: 1.0}
+
+	var final_result Results
+	final_result.New(n_images)
+
+	// Score each image by calculating the average score from camera burst
+	for i := range images {
+
+		//Evaluate camera burst
+		sample_result, err := ts.Evaluate(n_samples, images[i])
+
+		if err != nil {
+			fmt.Println(err)
+		}
+
+		// Append the average result from camera burst
+		final_result.Append(Average(sample_result.saturation), Average(sample_result.contrast), float64((i+1)*10), i)
+	}
+
+	// Plot the final results
+	err := PlotResults(final_result.turbidity, Normalize(final_result.saturation), Normalize(final_result.contrast))
+
+	fmt.Println(final_result.saturation)
+	fmt.Println(final_result.contrast)
 
-	err := ts.Evaluate(n, images)
 	if err != nil {
 		fmt.Println(err)
 	}
 }
 
-// RESULTS
-type Results struct {
-	turbidity  []float64
-	saturation []float64
-	contrast   []float64
-}
+// PLOTTING FUNCTIONS
 
-func (r *Results) New(n int) {
-	r.turbidity = make([]float64, n)
-	r.saturation = make([]float64, n)
-	r.contrast = make([]float64, n)
-}
-
-func (r *Results) Append(saturation, contrast, turbidity float64, i int) {
-	r.saturation[i] = saturation
-	r.contrast[i] = contrast
-	r.turbidity[i] = turbidity
-}
-
-func Normalize(slice []float64, n int) []float64 {
+// Normalize values in a slice between 0 and 1
+func Normalize(slice []float64) []float64 {
 
 	max := -1e10
 	min := 1e10
 
-	out := make([]float64, n)
+	out := make([]float64, len(slice))
 
-	if n <= 1 {
+	if len(slice) <= 1 {
 		return slice
 	}
 
-	for i := 0; i < n; i++ {
+	for i := range slice {
 		if slice[i] > max {
 			max = slice[i]
 		}
@@ -72,19 +83,30 @@ func Normalize(slice []float64, n int) []float64 {
 		}
 	}
 
-	for i := 0; i < n; i++ {
+	for i := range slice {
 		out[i] = (slice[i] - min) / (max - min)
 	}
 
 	return out
 }
 
-// PLOTTING
+// Return the average of a slice
+func Average(slice []float64) float64 {
+
+	out := 0.0
+
+	for i := range slice {
+		out += slice[i]
+	}
+
+	return out / float64(len(slice))
+}
 
 func PlotResults(x, saturation, contrast []float64) error {
+
 	err := plotToFile(
 		"Results",
-		"Filter Size",
+		"Turbidity (Almond Milk) (ml)",
 		"Score",
 		func(p *plot.Plot) error {
 			return plotutil.AddLinePoints(p,
@@ -93,9 +115,11 @@ func PlotResults(x, saturation, contrast []float64) error {
 			)
 		},
 	)
+
 	if err != nil {
 		return fmt.Errorf("Could not plot results: %v", err)
 	}
+
 	return nil
 }
 
@@ -121,7 +145,9 @@ func plotToFile(name, xTitle, yTitle string, draw func(*plot.Plot) error) error
 
 // plotterXY provides a plotter.XYs type value based on the given x and y data.
 func plotterXY(x, y []float64) plotter.XYs {
+
 	xy := make(plotter.XYs, len(x))
+
 	for i := range x {
 		xy[i].X = x[i]
 		xy[i].Y = y[i]
diff --git a/turbidity/results.go b/turbidity/results.go
index be3f58b2..701e037a 100644
--- a/turbidity/results.go
+++ b/turbidity/results.go
@@ -12,8 +12,8 @@ func (r *Results) New(n int) {
 	r.contrast = make([]float64, n)
 }
 
-func (r *Results) Append(saturation, contrast, turbidity float64, i int) {
-	r.saturation[i] = saturation
-	r.contrast[i] = contrast
-	r.turbidity[i] = turbidity
-}
\ No newline at end of file
+func (r *Results) Append(saturation, contrast, turbidity float64, index int) {
+	r.saturation[index] = saturation
+	r.contrast[index] = contrast
+	r.turbidity[index] = turbidity
+}
diff --git a/turbidity/turbidity.go b/turbidity/turbidity.go
index f7e754fa..cba8ebfe 100644
--- a/turbidity/turbidity.go
+++ b/turbidity/turbidity.go
@@ -2,20 +2,179 @@ package main
 
 import (
 	"fmt"
+	"image"
+	"math"
 
 	"gocv.io/x/gocv"
 )
 
-func main() {
-	var flag gocv.IMReadFlag = 0
-
-	window := gocv.NewWindow("Test")
-	img := gocv.IMRead("secci.jpg", flag)
-
-	fmt.Print(img.Size())
-
-	for {
-		window.IMShow(img)
-		window.WaitKey(0)
-	}
+// TURBIDITY SENSOR
+type TurbiditySensor struct {
+	template, standard gocv.Mat
+	k1, k2             int //block size
+	alpha, scale       float64
+}
+
+// Given a slice of test images of size n, return the sharpness and contrast scores
+func (ts TurbiditySensor) Evaluate(n int, images []gocv.Mat) (Results, error) {
+
+	var result Results
+	result.New(n)
+
+	for i := range images {
+
+		// Transform image
+		marker, err := ts.Transform(images[i])
+
+		if err != nil {
+			return result, fmt.Errorf("Image %v: %v", i, err)
+		}
+
+		// Apply sobel filter
+		edge := ts.Sobel(marker, ts.scale)
+
+		// Evaluate
+		out, err := ts.EvaluateImage(marker, edge, ts.k1, ts.k2, ts.alpha)
+
+		if err != nil {
+			return result, err
+		}
+
+		result.Append(out[0], out[1], float64((i+1)*10), i)
+	}
+
+	return result, nil
+}
+
+// Evaluate image sharpness and contrast using blocks of size k1 by k2. Return a slice of the respective scores
+func (ts TurbiditySensor) EvaluateImage(img, edge gocv.Mat, k1, k2 int, alpha float64) ([]float64, error) {
+
+	// Slice to store results
+	result := make([]float64, 2) // [0.0, 0.0]
+
+	if img.Rows()%k1 != 0 || img.Cols()%k2 != 0 {
+		return result, fmt.Errorf("Dimensions not compatible (%v, %v)", k1, k2)
+	}
+
+	l_step := int(img.Rows() / k1)
+	k_step := int(img.Cols() / k2)
+
+	for l := 0; l < img.Rows(); l = l + l_step {
+		for k := 0; k < img.Cols(); k = k + k_step {
+
+			//EME
+			err := ts.EvaluateBlock(edge, l, k, l+l_step, k+k_step, result, "EME", alpha)
+
+			if err != nil {
+				return result, err
+			}
+
+			//AMEE
+			err = ts.EvaluateBlock(img, l, k, l+l_step, k+k_step, result, "AMEE", alpha)
+
+			if err != nil {
+				return result, err
+			}
+		}
+	}
+
+	// EME
+	result[0] = 2.0 / (float64(k1) * float64(k1)) * result[0]
+
+	// AMEE
+	result[1] = -1.0 / (float64(k1) * float64(k1)) * result[1]
+
+	return result, nil
+}
+
+// Evaluate a block within and image and add to to the result slice
+func (ts TurbiditySensor) EvaluateBlock(img gocv.Mat, l1, k1, l2, k2 int, result []float64, operation string, alpha float64) error {
+
+	max := -1e10
+	min := 1e10
+
+	for i := l1; i < l2; i++ {
+		for j := k1; j < k2; j++ {
+
+			value := float64(img.GetUCharAt(i, j))
+
+			// Check max/min conditions, zero values are ignored
+			if value > max && value != 0.0 {
+				max = value
+			}
+			if value < min && value != 0.0 {
+				min = value
+			}
+		}
+	}
+
+	// Blocks which have no information are ignored
+	if max != -1e10 && min != 1e10 && max != min {
+		if operation == "EME" {
+
+			result[0] += math.Log(max / min)
+
+		} else if operation == "AMEE" {
+
+			contrast := (max + min) / (max - min)
+			result[1] += math.Pow(alpha*(contrast), alpha) * math.Log(contrast)
+
+		} else {
+			return fmt.Errorf("Invalid operation: %v", operation)
+		}
+	}
+	return nil
+}
+
+// Search image for matching template. Returns the transformed image  which best match the template
+func (ts TurbiditySensor) Transform(img gocv.Mat) (gocv.Mat, error) {
+
+	out := gocv.NewMat()
+	mask := gocv.NewMat()
+
+	corners_img := gocv.NewMat()
+	corners_template := gocv.NewMat()
+
+	// Find corners in image
+	if !gocv.FindChessboardCorners(img, image.Pt(3, 3), &corners_img, gocv.CalibCBNormalizeImage) {
+
+		// Apply default if transformation fails
+		// fmt.Println("Corner detection failed applying standard transformation")
+
+		if !gocv.FindChessboardCorners(ts.standard, image.Pt(3, 3), &corners_img, gocv.CalibCBNormalizeImage) {
+			return out, fmt.Errorf("Could not find corners in default image")
+		}
+	}
+
+	// Find corners in template
+	if !gocv.FindChessboardCorners(ts.template, image.Pt(3, 3), &corners_template, gocv.CalibCBNormalizeImage) {
+		return out, fmt.Errorf("Could not find corners in template")
+	}
+
+	// Find and apply transformation
+	H := gocv.FindHomography(corners_img, &corners_template, gocv.HomograpyMethodRANSAC, 3.0, &mask, 2000, 0.995)
+	gocv.WarpPerspective(img, &out, H, image.Pt(ts.template.Rows(), ts.template.Cols()))
+
+	return out, nil
+}
+
+// Apply sobel filter to an image with a given scale and return the result
+func (ts TurbiditySensor) Sobel(img gocv.Mat, scale float64) gocv.Mat {
+
+	dx := gocv.NewMat()
+	dy := gocv.NewMat()
+	sobel := gocv.NewMat()
+
+	// Apply filter
+	gocv.Sobel(img, &dx, gocv.MatTypeCV64F, 0, 1, 3, scale, 0.0, gocv.BorderConstant)
+	gocv.Sobel(img, &dy, gocv.MatTypeCV64F, 1, 0, 3, scale, 0.0, gocv.BorderConstant)
+
+	// Convert to unsigned
+	gocv.ConvertScaleAbs(dx, &dx, 1.0, 0.0)
+	gocv.ConvertScaleAbs(dy, &dy, 1.0, 0.0)
+
+	// Add x and y components
+	gocv.AddWeighted(dx, 0.5, dy, 0.5, 0, &sobel)
+
+	return sobel
 }