speed up SpatialMap smooth() in 2D without periodicity

Author: bhaller

VERSIONS

@@ -6,6 +6,8 @@ Note that not every commit will be logged here; that is what the Github commit h
 
 
 development head (in the master branch):
+	add cppcheck support in our cmake config, for development testing purposes
+	speed up the SpatialMap smooth() method in 2D without periodicity, since this is important for the Chevy et al. competition algorithm
 
 
 version 5.1 (Eidos version 4.1):

core/spatial_map.cpp

@@ -879,68 +879,184 @@ void SpatialMap::Convolve_S2(SpatialKernel &kernel)
 	double *new_values_ptr = new_values;
 
 	// FIXME: TO BE PARALLELIZED
-	for (int64_t b = 0; b < dim_b; ++b)
+	if (!periodic_a_ && !periodic_b_)
 	{
-		double coverage_b = (!periodic_b_ && ((b == 0) || (b == dim_b - 1))) ? 0.5 : 1.0;
-		
-		for (int64_t a = 0; a < dim_a; ++a)
+		// special-case for the completely non-periodic case with a non-zero kernel total
+		// see below for code that checks the optimizations done here, since they're a bit tricky
+		// FIXME this optimization should be extended to the 1D and 3D cases
+		for (int64_t b = 0; b < dim_b; ++b)
 		{
-			double coverage_a = (!periodic_a_ && ((a == 0) || (a == dim_a - 1))) ? 0.5 : 1.0;
-			double coverage = coverage_a * coverage_b;					// handles partial coverage at the edges of the spatial map
+			double coverage_b = ((b == 0) || (b == dim_b - 1)) ? 0.5 : 1.0;
 
-			// calculate the kernel's effect at point (a,b)
-			double kernel_total = 0.0;
-			double conv_total = 0.0;
+			// we want to loop over just the kernel extent that falls inside [0, dim_b - 1]
+			int64_t kernel_b_first = std::max((int64_t)0, -(b + kernel_b_offset));
+			int64_t kernel_b_last = std::min(kernel_dim_b - 1, (kernel_dim_b - 1) + (((dim_b - 1) - b) + kernel_b_offset));
 
-			for (int64_t kernel_a = 0; kernel_a < kernel_dim_a; kernel_a++)
+			for (int64_t a = 0; a < dim_a; ++a)
 			{
-				int64_t conv_a = a + kernel_a + kernel_a_offset;
+				double coverage_a = ((a == 0) || (a == dim_a - 1)) ? 0.5 : 1.0;
+				double coverage = coverage_a * coverage_b;					// handles partial coverage at the edges of the spatial map
+				
+				// we want to loop over just the kernel extent that falls inside [0, dim_a - 1]
+				int64_t kernel_a_first = std::max((int64_t)0, 0 - ((a - 0) + kernel_a_offset));
+				int64_t kernel_a_last = std::min(kernel_dim_a - 1, (kernel_dim_a - 1) + (((dim_a - 1) - a) + kernel_a_offset));
+				
+				// calculate the kernel's effect at point (a,b)
+				// BCH 10/2/2025: I tried adding a special case here for the middle portion where the kernel
+				// is not clipped at all and the kernel_total value is constant; it was not worthwhile,
+				// testing for whether we are in that case took time, and I think totalling the kernel
+				// in the loop below pretty much pipelines out so getting rid of it is inconsequential.
+				double kernel_total = 0.0;
+				double conv_total = 0.0;
 
-				// handle bounds: either clip or wrap
-				if ((conv_a < 0) || (conv_a >= dim_a))
+				for (int64_t kernel_a = kernel_a_first; kernel_a <= kernel_a_last; kernel_a++)
 				{
-					if (!periodic_a_)
-						continue;
+					int64_t conv_a = a + kernel_a + kernel_a_offset;
 
-					// periodicity: assume the two edges have identical values, skip over the edge value on the opposite side
-					while (conv_a < 0)
-						conv_a += (dim_a - 1);	// move -1 to dim - 2
-					while (conv_a >= dim_a)
-						conv_a -= (dim_a - 1);	// move dim to 1
+					for (int64_t kernel_b = kernel_b_first; kernel_b <= kernel_b_last; kernel_b++)
+					{
+						int64_t conv_b = b + kernel_b + kernel_b_offset;
+						
+						// this point is within bounds; add it in to the totals
+						double kernel_value = kernel_values[kernel_a + kernel_b * kernel_dim_a] * coverage;
+						double pixel_value = values_[conv_a + conv_b * dim_a];
+						
+						// we keep a total of the kernel values that were within bounds, for this point
+						kernel_total += kernel_value;
+						
+						// and we keep a total of the convolution - kernel values times pixel values
+						conv_total += kernel_value * pixel_value;
+					}
 				}
 
-				for (int64_t kernel_b = 0; kernel_b < kernel_dim_b; kernel_b++)
+				*(new_values_ptr++) = ((kernel_total > 0) ? (conv_total / kernel_total) : 0);
+			}
+		}
+	}
+	else
+	{
+		// general case that handles periodicity
+		for (int64_t b = 0; b < dim_b; ++b)
+		{
+			double coverage_b = (!periodic_b_ && ((b == 0) || (b == dim_b - 1))) ? 0.5 : 1.0;
+			
+			for (int64_t a = 0; a < dim_a; ++a)
+			{
+				double coverage_a = (!periodic_a_ && ((a == 0) || (a == dim_a - 1))) ? 0.5 : 1.0;
+				double coverage = coverage_a * coverage_b;					// handles partial coverage at the edges of the spatial map
+				
+				// calculate the kernel's effect at point (a,b)
+				double kernel_total = 0.0;
+				double conv_total = 0.0;
+				
+				for (int64_t kernel_a = 0; kernel_a < kernel_dim_a; kernel_a++)
 				{
-					int64_t conv_b = b + kernel_b + kernel_b_offset;
+					int64_t conv_a = a + kernel_a + kernel_a_offset;
 
 					// handle bounds: either clip or wrap
-					if ((conv_b < 0) || (conv_b >= dim_b))
+					if ((conv_a < 0) || (conv_a >= dim_a))
 					{
-						if (!periodic_b_)
+						if (!periodic_a_)
 							continue;
 
 						// periodicity: assume the two edges have identical values, skip over the edge value on the opposite side
-						while (conv_b < 0)
-							conv_b += (dim_b - 1);	// move -1 to dim - 2
-						while (conv_b >= dim_b)
-							conv_b -= (dim_b - 1);	// move dim to 1
+						while (conv_a < 0)
+							conv_a += (dim_a - 1);	// move -1 to dim - 2
+						while (conv_a >= dim_a)
+							conv_a -= (dim_a - 1);	// move dim to 1
 					}
 
-					// this point is within bounds; add it in to the totals
-					double kernel_value = kernel_values[kernel_a + kernel_b * kernel_dim_a] * coverage;
-					double pixel_value = values_[conv_a + conv_b * dim_a];
+					for (int64_t kernel_b = 0; kernel_b < kernel_dim_b; kernel_b++)
+					{
+						int64_t conv_b = b + kernel_b + kernel_b_offset;
+						
+						// handle bounds: either clip or wrap
+						if ((conv_b < 0) || (conv_b >= dim_b))
+						{
+							if (!periodic_b_)
+								continue;
+							
+							// periodicity: assume the two edges have identical values, skip over the edge value on the opposite side
+							while (conv_b < 0)
+								conv_b += (dim_b - 1);	// move -1 to dim - 2
+							while (conv_b >= dim_b)
+								conv_b -= (dim_b - 1);	// move dim to 1
+						}
+						
+						// this point is within bounds; add it in to the totals
+						double kernel_value = kernel_values[kernel_a + kernel_b * kernel_dim_a] * coverage;
+						double pixel_value = values_[conv_a + conv_b * dim_a];
+						
+						// we keep a total of the kernel values that were within bounds, for this point
+						kernel_total += kernel_value;
+						
+						// and we keep a total of the convolution - kernel values times pixel values
+						conv_total += kernel_value * pixel_value;
+					}
+				}
+				
+				*(new_values_ptr++) = ((kernel_total > 0) ? (conv_total / kernel_total) : 0);
+			}
+		}
+	}
+	
+#if 0
+	// Check the values computed above by the optimized algorithm against the simple algorithm
+	// test_smooth2D.slim is a test file that can be used with this validation code
+	new_values_ptr = new_values;
+	
+	if (!periodic_a_ && !periodic_b_)
+	{
+		for (int64_t b = 0; b < dim_b; ++b)
+		{
+			double coverage_b = (!periodic_b_ && ((b == 0) || (b == dim_b - 1))) ? 0.5 : 1.0;
+			
+			for (int64_t a = 0; a < dim_a; ++a)
+			{
+				double coverage_a = (!periodic_a_ && ((a == 0) || (a == dim_a - 1))) ? 0.5 : 1.0;
+				double coverage = coverage_a * coverage_b;					// handles partial coverage at the edges of the spatial map
+				
+				// calculate the kernel's effect at point (a,b)
+				double kernel_total = 0.0;
+				double conv_total = 0.0;
+				
+				for (int64_t kernel_a = 0; kernel_a < kernel_dim_a; kernel_a++)
+				{
+					int64_t conv_a = a + kernel_a + kernel_a_offset;
 
-					// we keep a total of the kernel values that were within bounds, for this point
-					kernel_total += kernel_value;
+					// handle bounds: either clip or wrap
+					if ((conv_a < 0) || (conv_a >= dim_a))
+						continue;
 
-					// and we keep a total of the convolution - kernel values times pixel values
-					conv_total += kernel_value * pixel_value;
+					for (int64_t kernel_b = 0; kernel_b < kernel_dim_b; kernel_b++)
+					{
+						int64_t conv_b = b + kernel_b + kernel_b_offset;
+						
+						// handle bounds: either clip or wrap
+						if ((conv_b < 0) || (conv_b >= dim_b))
+							continue;
+						
+						// this point is within bounds; add it in to the totals
+						double kernel_value = kernel_values[kernel_a + kernel_b * kernel_dim_a] * coverage;
+						double pixel_value = values_[conv_a + conv_b * dim_a];
+						
+						// we keep a total of the kernel values that were within bounds, for this point
+						kernel_total += kernel_value;
+						
+						// and we keep a total of the convolution - kernel values times pixel values
+						conv_total += kernel_value * pixel_value;
+					}
 				}
+				
+				double check_value = ((kernel_total > 0) ? (conv_total / kernel_total) : 0);
+				double calculated_value = *(new_values_ptr++);
+				
+				if (check_value != calculated_value)
+					std::cout << "Convolve_S2 optimization mismatch: " << check_value << " versus " << calculated_value << std::endl;
 			}
-			
-			*(new_values_ptr++) = ((kernel_total > 0) ? (conv_total / kernel_total) : 0);
 		}
 	}
+#endif
 
 	TakeOverMallocedValues(new_values, 2, grid_size_);	// takes new_values from us
 }