NLM searchWindow = xy.z for anisotropy

src/MATLAB/+HIP/@Cuda/Cuda.m

@@ -7,13 +7,13 @@ methods (Static)
     [hydraConfig] = CheckConfig()
     [imageOut] = IdentityFilter(imageIn,device)
     [cmdInfo] = Info()
-    [imageOut] = Gaussian(imageIn,sigmas,numIterations,device)
-    [imageOut] = MeanFilter(imageIn,kernel,numIterations,device)
     [deviceStatsArray] = DeviceStats()
     [imageOut] = ElementWiseDifference(image1In,image2In,device)
     [imageOut] = MultiplySum(imageIn,kernel,numIterations,device)
     [imageOut] = EntropyFilter(imageIn,kernel,device)
     [imageOut] = StdFilter(imageIn,kernel,numIterations,device)
+    [imageOut] = Gaussian(imageIn,sigmas,numIterations,device)
+    [imageOut] = MeanFilter(imageIn,kernel,numIterations,device)
     [minVal,maxVal] = GetMinMax(imageIn,device)
     [imageOut] = LoG(imageIn,sigmas,device)
     [imageOut] = MedianFilter(imageIn,kernel,numIterations,device)

src/c/Cuda/CudaNLMeans.cuh

@@ -19,29 +19,32 @@
 template <class PixelTypeIn, class PixelTypeOut>
 // params needed: a,h, search window size, comparison nhood
 __global__ void cudaNLMeans_mv(CudaImageContainer<PixelTypeIn> imageIn, CudaImageContainer<PixelTypeOut> imageMean, CudaImageContainer<PixelTypeOut> imageVariance,
-	CudaImageContainer<PixelTypeOut> imageOut, double h, int searchWindowRadius, int nhoodRadius, PixelTypeOut minValue, PixelTypeOut maxValue)
+	CudaImageContainer<PixelTypeOut> imageOut, double h, double searchWindowRadius, int nhoodRadius, PixelTypeOut minValue, PixelTypeOut maxValue)
 {
 	Vec<std::size_t> threadCoordinate;
 	GetThreadBlockCoordinate(threadCoordinate);
 
 	if (threadCoordinate < imageIn.getDims())
 	{
-		Vec<int> searchVec = Vec<int>::min(Vec<int>(imageIn.getDims()) - 1, Vec<int>(searchWindowRadius));
+		int xySearch = (int)std::floor(searchWindowRadius);
+		int zSearch = (int)((searchWindowRadius - std::floor(searchWindowRadius)) * (double)xySearch);
+		if ( (searchWindowRadius - std::floor(searchWindowRadius)) < 1e-6 )
+			zSearch = xySearch;
+		Vec<int> nSearchWindow = Vec<int>(xySearch, xySearch, zSearch);
+		Vec<int> searchVec = Vec<int>::min(Vec<int>(imageIn.getDims()) - 1, nSearchWindow);
 
 		// 
 		float inputVal = (float)imageIn(threadCoordinate);
 		float inputMeanVal = (float)imageMean(threadCoordinate);
 		float inputVarVal = (float)imageVariance(threadCoordinate);
 		// 
-		double wMax = 0.;
+		double wMax = 0.; 
 		double wAccumulator = 0.;
 		double outputAccumulator = 0.;
 		KernelIterator kIt(threadCoordinate, imageIn.getDims(), searchVec * 2 + 1);
 		for (; !kIt.end(); ++kIt)
 		{
 			Vec<float> kernelPos = kIt.getImageCoordinate();
-			if (kernelPos == threadCoordinate)
-				continue;
 
 			float kernelMeanVal = (float)imageMean(kernelPos);
 			float kernelVarVal = (float)imageVariance(kernelPos);
@@ -65,84 +68,11 @@ __global__ void cudaNLMeans_mv(CudaImageContainer<PixelTypeIn> imageIn, CudaImag
 	}
 } // cudaNLMeans_mv
 
-template <class PixelTypeIn, class PixelTypeOut>
-// params needed: a,h, search window size, comparison nhood
-__global__ void cudaNLMeans(CudaImageContainer<PixelTypeIn> imageIn, CudaImageContainer<PixelTypeOut> imageOut, double h, 
-	int searchWindowRadius, int nhoodRadius, PixelTypeOut minValue, PixelTypeOut maxValue)
-{
-	Vec<std::size_t> threadCoordinate;
-	GetThreadBlockCoordinate(threadCoordinate);
-
-	if (threadCoordinate < imageIn.getDims())
-	{
-		// for loop over whole image, except for searchWindowSize
-		Vec<int> itc = Vec<int>(threadCoordinate);
-
-		Vec<int> nhoodVec = Vec<int>::min(Vec<int>(imageIn.getDims())-1, Vec<int>(nhoodRadius));
-
-		Vec<int> searchMin = Vec<int>::max(itc - searchWindowRadius, nhoodVec);
-		Vec<int> searchMax = Vec<int>::min(itc + searchWindowRadius+1, Vec<int>(imageIn.getDims()) - nhoodVec);
-
-		// outValAccumulator is unnormalized output value -- we normalize at the end
-		double outValAccumulator= 0.0; // 
-		// wAccumulator is for the normalizing constant
-		double wAccumulator = 1e-7; // 
-		double wMax = 0; // for application when x==y==z==0
-		for (int z = searchMin.z; z < searchMax.z; z++) {
-			for (int y = searchMin.y; y < searchMax.y; y++) {
-				for (int x = searchMin.x; x < searchMax.x; x++) {
-										
-					if ((0 == x) && (0 == y) && (0 == z))
-						continue;
-
-					// center the kernel on threadCoordinate
-					// iterator over nhoodSize
-					KernelIterator kIt(threadCoordinate, imageIn.getDims(), nhoodVec*2+1);
-					// image value at the x,y,z neighborhood window center
-					float searchVal = (float)imageIn( Vec<float>(x,y,z));
-					float wi=0.;
-					for (; !kIt.end(); ++kIt)
-					{
-						Vec<float> imInPos = kIt.getImageCoordinate();
-						float inVal = (float)imageIn(imInPos);
-						 
-						// now get the corresponding pixel value in the search window space
-						Vec<float> sCoord = Vec<float>(x,y,z)-nhoodVec+ kIt.getKernelCoordinate();
-						float sVal = (float)imageIn(sCoord);
-
-						double diff = SQR(inVal - sVal);
-						// weight by gaussian (Ga from paper)
-						double ga = 1; // exp(-1 * (Vec<int>(kIt.getKernelCoordinate()) - nhoodVec).lengthSqr() / SQR(a));
-						// w is gaussian weighted euclidean distance
-						wi += ga*diff;
-					}
-					double w = exp(-wi / SQR(h));
-
-					if (w > wMax)
-						wMax = w;
-
-					outValAccumulator += w*searchVal;
-					wAccumulator += w;
-				}
-			}
-		}
-		// add in for x==y==z==0
-		outValAccumulator += wMax*(float)imageIn(threadCoordinate);
-		wAccumulator += wMax;
-
-
-		// now normalize
-		double outVal = outValAccumulator / wAccumulator;
-				
-		imageOut(threadCoordinate) = (PixelTypeOut)CLAMP(outVal, minValue, maxValue);
-	}
-} // cudaNLMeans 
-
 // this part is the templated function that gets called by the front end.
 // here be cpu
 // todo - if we chunk, make sure search window doesn't go off chunk
 template <class PixelTypeIn, class PixelTypeOut>
-void cNLMeans(ImageView<PixelTypeIn> imageIn, ImageView<PixelTypeOut> imageOut, double h, int searchWindowRadius, int nhoodRadius, int device = -1)
+void cNLMeans(ImageView<PixelTypeIn> imageIn, ImageView<PixelTypeOut> imageOut, double h, double searchWindowRadius, int nhoodRadius, int device = -1)
 {
 	const PixelTypeOut MIN_VAL = std::numeric_limits<PixelTypeOut>::lowest();
 	const PixelTypeOut MAX_VAL = std::numeric_limits<PixelTypeOut>::max();