diff --git a/src/dsp/gpu/ConvGpu.java b/src/dsp/gpu/ConvGpu.java index ffd0af2..cd3f10a 100644 --- a/src/dsp/gpu/ConvGpu.java +++ b/src/dsp/gpu/ConvGpu.java @@ -32,7 +32,7 @@ public class ConvGpu implements IConv { JCudaDriver.setExceptionsEnabled(true); JCudaDriver.cuInit(0); } - + public ConvGpu() { try { // Initialize CUDA context @@ -70,8 +70,60 @@ public ConvGpu() { } } + private float[] convolve2DGpu(float[] input, int width, int height, float[] kernel2D, int kernelSizeInt) { + float[] output = new float[input.length]; + + long inputSizeBytes = (long) input.length * Sizeof.FLOAT; + long kernelSizeBytes = (long) kernel2D.length * Sizeof.FLOAT; + long outputSizeBytes = (long) output.length * Sizeof.FLOAT; + + CUdeviceptr d_input = null; + CUdeviceptr d_kernel = null; + CUdeviceptr d_output = null; + + try { + d_input = CudaMemoryManager.allocate(inputSizeBytes); + d_kernel = CudaMemoryManager.allocate(kernelSizeBytes); + d_output = CudaMemoryManager.allocate(outputSizeBytes); + + JCudaDriver.cuMemcpyHtoDAsync(d_input, Pointer.to(input), inputSizeBytes, stream); + JCudaDriver.cuMemcpyHtoDAsync(d_kernel, Pointer.to(kernel2D), kernelSizeBytes, stream); + + Pointer kernelParameters = Pointer.to( + Pointer.to(d_input), + Pointer.to(d_kernel), + Pointer.to(d_output), + Pointer.to(new int[]{width}), + Pointer.to(new int[]{height}), + Pointer.to(new int[]{kernelSizeInt}) + ); + + int blockSize = 16; + int gridX = (width + blockSize - 1) / blockSize; + int gridY = (height + blockSize - 1) / blockSize; + + JCudaDriver.cuLaunchKernel( + function2D, + gridX, gridY, 1, + blockSize, blockSize, 1, + 0, stream, + kernelParameters, + null + ); + + JCudaDriver.cuMemcpyDtoHAsync(Pointer.to(output), d_output, outputSizeBytes, stream); + JCudaDriver.cuStreamSynchronize(stream); + + return output; + } finally { + if (d_input != null) CudaMemoryManager.free(d_input, inputSizeBytes); + if (d_kernel != null) CudaMemoryManager.free(d_kernel, kernelSizeBytes); + if (d_output != null) CudaMemoryManager.free(d_output, outputSizeBytes); + } + } + @Override - public void convolveSemiSep(FloatProcessor ip, float[] kernx, float[] kern_diff) { + public void convolveSemiSep(FloatProcessor ip, float[] kernx, float[] kern_diff) { FloatProcessor ip2 = null; FloatProcessor ipx = null; final Rectangle roi = ip.getRoi(); @@ -85,19 +137,20 @@ public void convolveSemiSep(FloatProcessor ip, float[] kernx, float[] kern_diff) ipx.setSnapshotPixels(ip.getSnapshotPixels()); } - convolveFloat1D(ipx, kern_diff, kern_diff.length, 1); // x direction + convolveFloat1D(ipx, kern_diff, Ox); ipx.setSnapshotPixels(null); - convolveFloat1D(ipx, kernx, 1, kernx.length); // y direction + convolveFloat1D(ipx, kernx, Oy); - convolveFloat1D(ip2, kernx, kernx.length, 1); // x direction + convolveFloat1D(ip2, kernx, Ox); ip2.setSnapshotPixels(null); - convolveFloat1D(ip2, kern_diff, 1, kern_diff.length); // y direction + convolveFloat1D(ip2, kern_diff, Oy); + add(ip2, ipx, ip2.getRoi()); ip.setPixels(ip2.getPixels()); } @Override - public void convolveSemiSepIter(FloatProcessor ip, float[] kernx, float[] kern_diff) { + public void convolveSemiSepIter(FloatProcessor ip, float[] kernx, float[] kern_diff) { FloatProcessor ip2 = (FloatProcessor) ip.duplicate(); FloatProcessor ipx = (FloatProcessor) ip.duplicate(); final Rectangle roi = ip.getRoi(); @@ -105,30 +158,32 @@ public void convolveSemiSepIter(FloatProcessor ip, float[] kernx, float[] kern_d ip2.setRoi(roi); ipx.setRoi(roi); - convolveFloat1D(ipx, kern_diff, Ox); // x direction - convolveFloat1D(ipx, kernx, Oy); // y direction + convolveFloat1D(ipx, kern_diff, Ox); + convolveFloat1D(ipx, kernx, Oy); - convolveFloat1D(ip2, kernx, Ox); // x direction - convolveFloat1D(ip2, kern_diff, Oy); // y direction + convolveFloat1D(ip2, kernx, Ox); + convolveFloat1D(ip2, kern_diff, Oy); add(ip2, ipx, ip.getRoi()); ip.setPixels(ip2.getPixels()); } @Override - public void convolveSepIter(FloatProcessor ip, float[] kernx, float[] kern_diff) { - convolveFloat1D(ip, kern_diff, Ox); // x direction - convolveFloat1D(ip, kernx, Oy); // y direction + public void convolveSepIter(FloatProcessor ip, float[] kernx, float[] kern_diff) { + convolveFloat1D(ip, kern_diff, Ox); + convolveFloat1D(ip, kernx, Oy); } @Override - public void convolveSep(ImageProcessor ip, float[] kernx, float[] kern_diff) { - convolveFloat1D(ip, kern_diff, kern_diff.length, 1); // x direction - convolveFloat1D(ip, kernx, 1, kernx.length); // y direction + public void convolveSep(ImageProcessor ip, float[] kernx, float[] kern_diff) { + // FIX: convolveSep must be separable 1D X then 1D Y. + // It should NOT call the kw/kh overload that uses the 2D CUDA kernel. + convolveFloat1D((FloatProcessor) ip, kern_diff, Ox); // X + convolveFloat1D((FloatProcessor) ip, kernx, Oy); // Y } @Override - public void convolveSemiSep(ImageStack xstack, float[] kernx, float[] kerny, float[] kernz) { + public void convolveSemiSep(ImageStack xstack, float[] kernx, float[] kerny, float[] kernz) { long time = -System.nanoTime(); ImageStack ystack = cloneStack(xstack); ImageStack zstack = cloneStack(xstack); @@ -137,17 +192,17 @@ public void convolveSemiSep(ImageStack xstack, float[] kernx, float[] kerny, flo time /= 1000.0f; time = -System.nanoTime(); - convolveFloat1D(xstack, kernx, Ox); // X - convolveFloat1D(xstack, kerny, Oy); // Y - convolveFloat1D(xstack, kernz, Oz); // Z + convolveFloat1D(xstack, kernx, Ox); + convolveFloat1D(xstack, kerny, Oy); + convolveFloat1D(xstack, kernz, Oz); - convolveFloat1D(ystack, kernx, Oy); // Y - convolveFloat1D(ystack, kerny, Ox); // X - convolveFloat1D(ystack, kernz, Oz); // Z + convolveFloat1D(ystack, kernx, Oy); + convolveFloat1D(ystack, kerny, Ox); + convolveFloat1D(ystack, kernz, Oz); - convolveFloat1D(zstack, kernx, Oz); // Z - convolveFloat1D(zstack, kerny, Ox); // X - convolveFloat1D(zstack, kernz, Oy); // Y + convolveFloat1D(zstack, kernx, Oz); + convolveFloat1D(zstack, kerny, Ox); + convolveFloat1D(zstack, kernz, Oy); addToStack(xstack, ystack, zstack); ystack = null; @@ -159,7 +214,7 @@ public void convolveSemiSep(ImageStack xstack, float[] kernx, float[] kerny, flo } @Override - public void convolveSep3D(ImageStack xstack, float[] kernx, float[] kern_diffx, float[] kernz) { + public void convolveSep3D(ImageStack xstack, float[] kernx, float[] kern_diffx, float[] kernz) { convolveFloat1D(xstack, kern_diffx, Ox); convolveFloat1D(xstack, kernx, Oy); convolveFloat1D(xstack, kernz, Oz); @@ -167,8 +222,7 @@ public void convolveSep3D(ImageStack xstack, float[] kernx, float[] kern_diffx, private void addToStack(ImageStack dest, ImageStack a, ImageStack b) { int bitdepth = dest.getBitDepth(); - if (bitdepth != a.getBitDepth() || a.getBitDepth() != b.getBitDepth()) - return; + if (bitdepth != a.getBitDepth() || a.getBitDepth() != b.getBitDepth()) return; final int sz = dest.getSize(); for (int i = 1; i <= sz; i++) { @@ -177,32 +231,28 @@ private void addToStack(ImageStack dest, ImageStack a, ImageStack b) { byte[] pixels = (byte[]) dest.getPixels(i); byte[] pixels_a = (byte[]) a.getPixels(i); byte[] pixels_b = (byte[]) b.getPixels(i); - for (int c = 0; c < pixels.length; c++) - pixels[c] += pixels_a[c] + pixels_b[c]; + for (int c = 0; c < pixels.length; c++) pixels[c] += pixels_a[c] + pixels_b[c]; break; } case 16: { short[] pixels = (short[]) dest.getPixels(i); short[] pixels_a = (short[]) a.getPixels(i); short[] pixels_b = (short[]) b.getPixels(i); - for (int c = 0; c < pixels.length; c++) - pixels[c] += pixels_a[c] + pixels_b[c]; + for (int c = 0; c < pixels.length; c++) pixels[c] += pixels_a[c] + pixels_b[c]; break; } case 24: { int[] pixels = (int[]) dest.getPixels(i); int[] pixels_a = (int[]) a.getPixels(i); int[] pixels_b = (int[]) b.getPixels(i); - for (int c = 0; c < pixels.length; c++) - pixels[c] += pixels_a[c] + pixels_b[c]; + for (int c = 0; c < pixels.length; c++) pixels[c] += pixels_a[c] + pixels_b[c]; break; } case 32: { float[] pixels = (float[]) dest.getPixels(i); float[] pixels_a = (float[]) a.getPixels(i); float[] pixels_b = (float[]) b.getPixels(i); - for (int c = 0; c < pixels.length; c++) - pixels[c] += pixels_a[c] + pixels_b[c]; + for (int c = 0; c < pixels.length; c++) pixels[c] += pixels_a[c] + pixels_b[c]; break; } } @@ -216,8 +266,7 @@ public static ImageStack cloneStack(ImageStack is) { ImageStack ret = ImageStack.create(width, height, array.length, is.getBitDepth()); Object[] array2 = array.clone(); int cnt = 1; - for (Object o : array2) - ret.setPixels(o, cnt++); + for (Object o : array2) ret.setPixels(o, cnt++); ret.update(is.getProcessor(1)); ret.setRoi(is.getRoi()); return ret; @@ -233,7 +282,7 @@ private void add(ImageProcessor dest, ImageProcessor src, Rectangle r) { } @Override - public boolean convolveFloat(ImageProcessor ip, float[] kernel, int kw, int kh) { + public boolean convolveFloat(ImageProcessor ip, float[] kernel, int kw, int kh) { if (!gpuInitialized) { IJ.log("GPU not initialized - cannot perform convolution"); return false; @@ -249,7 +298,6 @@ private boolean convolveFloatGPU(ImageProcessor ip, float[] kernel, int kw, int float[] output = new float[input.length]; try { - // Use memory pool long inputSize = input.length * Sizeof.FLOAT; long kernelSize = kernel.length * Sizeof.FLOAT; long outputSize = output.length * Sizeof.FLOAT; @@ -258,11 +306,9 @@ private boolean convolveFloatGPU(ImageProcessor ip, float[] kernel, int kw, int CUdeviceptr d_kernel = CudaMemoryManager.allocate(kernelSize); CUdeviceptr d_output = CudaMemoryManager.allocate(outputSize); - // Use async memory copies with stream JCudaDriver.cuMemcpyHtoDAsync(d_input, Pointer.to(input), inputSize, stream); JCudaDriver.cuMemcpyHtoDAsync(d_kernel, Pointer.to(kernel), kernelSize, stream); - // Setup kernel parameters Pointer kernelParameters = Pointer.to( Pointer.to(d_input), Pointer.to(d_kernel), @@ -272,7 +318,6 @@ private boolean convolveFloatGPU(ImageProcessor ip, float[] kernel, int kw, int Pointer.to(new int[]{kw}) ); - // Launch kernel with stream int blockSize = 16; int gridX = (width + blockSize - 1) / blockSize; int gridY = (height + blockSize - 1) / blockSize; @@ -282,16 +327,11 @@ private boolean convolveFloatGPU(ImageProcessor ip, float[] kernel, int kw, int blockSize, blockSize, 1, 0, stream, kernelParameters, null); - // Async copy back to host JCudaDriver.cuMemcpyDtoHAsync(Pointer.to(output), d_output, outputSize, stream); - - // Wait for all async operations to complete JCudaDriver.cuStreamSynchronize(stream); - // Set the result ip.setPixels(output); - // Return memory to pool CudaMemoryManager.free(d_input, inputSize); CudaMemoryManager.free(d_kernel, kernelSize); CudaMemoryManager.free(d_output, outputSize); @@ -308,10 +348,7 @@ private boolean convolveFloatGPU(ImageProcessor ip, float[] kernel, int kw, int public void cleanup() { try { - if (stream != null) { - JCudaDriver.cuStreamDestroy(stream); - } - // Free all pooled memory at application exit + if (stream != null) JCudaDriver.cuStreamDestroy(stream); CudaMemoryManager.freeAll(); } catch (Exception e) { IJ.log("GPU cleanup failed: " + e.getMessage()); @@ -319,34 +356,107 @@ public void cleanup() { } @Override - public void convolveFloat1D(FloatProcessor fp, float[] kernel, int xdir) { - IJLineIteratorIP iter = new IJLineIteratorIP(fp, xdir); + public void convolveFloat1D(FloatProcessor fp, float[] kernel, int xdir) { final int width = fp.getWidth(); final int height = fp.getHeight(); - FloatProcessor ret = new FloatProcessor(width, height); - int cnt = 0; - if (debug) { - printvector(kernel); - System.out.println(); + if (!gpuInitialized || kernel == null || kernel.length == 0 || (xdir != Ox && xdir != Oy)) { + // Fallback: original per-line implementation (still correct). + IJLineIteratorIP iter = new IJLineIteratorIP(fp, xdir); + FloatProcessor ret = new FloatProcessor(width, height); + int cnt = 0; + while (iter.hasNext()) { + final float[] line = iter.next(); + final float[] line2 = lineConvolveGPU(line, kernel, false); + iter.putLineFloat(ret, line2, cnt, xdir); + cnt++; + } + fp.setPixels(ret.getPixels()); + return; } - while (iter.hasNext()) { - final float[] line = iter.next(); - final float[] line2 = lineConvolveGPU(line, kernel, false); - iter.putLineFloat(ret, line2, cnt, xdir); - cnt++; + + final int kernelSizeInt = kernel.length; + final int halfKernel = kernelSizeInt / 2; + + // Build sparse "embedded" 2D kernel so convolve2DKernel behaves like a 1D convolution along axis. + float[] kernel2D = new float[kernelSizeInt * kernelSizeInt]; + if (xdir == Ox) { + int base = halfKernel * kernelSizeInt; + for (int col = 0; col < kernelSizeInt; col++) kernel2D[base + col] = kernel[col]; + } else { // Oy + int col = halfKernel; + for (int row = 0; row < kernelSizeInt; row++) kernel2D[row * kernelSizeInt + col] = kernel[row]; + } + + float[] paddedInput; + int paddedWidth; + int paddedHeight; + + if (xdir == Ox) { + paddedWidth = width + 2 * halfKernel; + paddedHeight = height; + paddedInput = new float[paddedWidth * paddedHeight]; + + float[] input = (float[]) fp.getPixels(); + for (int y = 0; y < height; y++) { + int srcRowOff = y * width; + int dstRowOff = y * paddedWidth; + + float leftVal = input[srcRowOff]; + for (int x = 0; x < halfKernel; x++) paddedInput[dstRowOff + x] = leftVal; + + System.arraycopy(input, srcRowOff, paddedInput, dstRowOff + halfKernel, width); + + float rightVal = input[srcRowOff + width - 1]; + for (int x = 0; x < halfKernel; x++) paddedInput[dstRowOff + halfKernel + width + x] = rightVal; + } + + float[] paddedOut = convolve2DGpu(paddedInput, paddedWidth, paddedHeight, kernel2D, kernelSizeInt); + + float[] out = new float[width * height]; + for (int y = 0; y < height; y++) { + int srcOff = y * paddedWidth + halfKernel; + int dstOff = y * width; + System.arraycopy(paddedOut, srcOff, out, dstOff, width); + } + fp.setPixels(out); + } else { // Oy + paddedWidth = width; + paddedHeight = height + 2 * halfKernel; + paddedInput = new float[paddedWidth * paddedHeight]; + + float[] input = (float[]) fp.getPixels(); + for (int py = 0; py < paddedHeight; py++) { + int sy = py - halfKernel; + if (sy < 0) sy = 0; + if (sy >= height) sy = height - 1; + + int srcRowOff = sy * width; + int dstRowOff = py * paddedWidth; + System.arraycopy(input, srcRowOff, paddedInput, dstRowOff, width); + } + + float[] paddedOut = convolve2DGpu(paddedInput, paddedWidth, paddedHeight, kernel2D, kernelSizeInt); + + float[] out = new float[width * height]; + for (int y = 0; y < height; y++) { + int srcOff = (y + halfKernel) * paddedWidth; + int dstOff = y * width; + System.arraycopy(paddedOut, srcOff, out, dstOff, width); + } + fp.setPixels(out); } - fp.setPixels(ret.getPixels()); } @Override - public void convolveFloat1D(ImageStack is, float[] kernel, int xdir) { + public void convolveFloat1D(ImageStack is, float[] kernel, int xdir) { IJLineIteratorStack iter = new IJLineIteratorStack(is, xdir); final int width = is.getWidth(); final int height = is.getHeight(); final int depth = is.getSize(); ImageStack ret = ImageStack.create(width, height, depth, is.getBitDepth()); int cnt = 0; + while (iter.hasNext()) { final float[] line = iter.next(); final float[] line2 = lineConvolveGPU(line, kernel, false); @@ -361,42 +471,33 @@ public void convolveFloat1D(ImageStack is, float[] kernel, int xdir) { } @Override - public void convolveFloat1D(ImageProcessor ip, float[] kernel, int kw, int kh) { + public void convolveFloat1D(ImageProcessor ip, float[] kernel, int kw, int kh) { convolveFloatGPU(ip, kernel, kw, kh); } static void printvector(float[] data) { - for (int i = 0; i < data.length; i++) { - System.out.print(data[i] + ","); - } + for (int i = 0; i < data.length; i++) System.out.print(data[i] + ","); } public float[] lineConvolveGPU(float[] arr, float[] kernel, boolean flip) { - if (!gpuInitialized || arr == null || kernel == null || arr.length == 0) { - return lineConvolve(arr, kernel, flip); - } + if (!gpuInitialized || arr == null || kernel == null || arr.length == 0) return lineConvolve(arr, kernel, flip); - long startTime = System.nanoTime(); int inputLength = arr.length; int kernelSize = kernel.length; float[] output = new float[inputLength]; try { - // Calculate memory sizes long inputSize = inputLength * Sizeof.FLOAT; long kernelSizeBytes = kernelSize * Sizeof.FLOAT; long outputSize = inputLength * Sizeof.FLOAT; - //allocate device memory using memory pool CUdeviceptr d_input = CudaMemoryManager.allocate(inputSize); CUdeviceptr d_kernel = CudaMemoryManager.allocate(kernelSizeBytes); CUdeviceptr d_output = CudaMemoryManager.allocate(outputSize); - // Copy data to device JCudaDriver.cuMemcpyHtoDAsync(d_input, Pointer.to(arr), inputSize, stream); JCudaDriver.cuMemcpyHtoDAsync(d_kernel, Pointer.to(kernel), kernelSizeBytes, stream); - // stup kernel parameters Pointer kernelParameters = Pointer.to( Pointer.to(d_input), Pointer.to(d_kernel), @@ -406,46 +507,32 @@ public float[] lineConvolveGPU(float[] arr, float[] kernel, boolean flip) { Pointer.to(new int[]{flip ? 1 : 0}) ); - //config launch param int blockSize = 256; int gridSize = (inputLength + blockSize - 1) / blockSize; - // Launch the kernel JCudaDriver.cuLaunchKernel(functionLineConvolve1D, - gridSize, 1, 1, // grid dim - blockSize, 1, 1, // block dim - 0, stream, // shared memory n stream - kernelParameters, null); + gridSize, 1, 1, + blockSize, 1, 1, + 0, stream, + kernelParameters, + null); - // Copy result back to host JCudaDriver.cuMemcpyDtoHAsync(Pointer.to(output), d_output, outputSize, stream); - - // Wait for all operations to complete JCudaDriver.cuStreamSynchronize(stream); - // Free memory back to pool CudaMemoryManager.free(d_input, inputSize); CudaMemoryManager.free(d_kernel, kernelSizeBytes); CudaMemoryManager.free(d_output, outputSize); - long endTime = System.nanoTime(); - if (debug) { - IJ.log(String.format("GPU 1D convolution time: %.2f ms (size: %d)", - (endTime - startTime) / 1e6, inputLength)); - } - return output; - } catch (Exception e) { IJ.log("GPU 1D convolution failed: " + e.getMessage()); - // Fall back to CPU implementation return lineConvolve(arr, kernel, flip); } } private static float[] lineConvolve(float[] arr, float[] kernel, boolean flip) { - if (flip) - flip(kernel); + if (flip) flip(kernel); float[] y = new float[arr.length]; int kw = kernel.length / 2; @@ -454,11 +541,8 @@ private static float[] lineConvolve(float[] arr, float[] kernel, boolean flip) { int c = 0; for (int k = -kw; k <= kw; k++) { int q = i - k; - if (0 <= q && q < arr.length) { - y[i] += arr[q] * kernel[c]; - } else { - y[i] += arr[0] * kernel[c]; - } + if (0 <= q && q < arr.length) y[i] += arr[q] * kernel[c]; + else y[i] += arr[0] * kernel[c]; c++; } } @@ -475,14 +559,12 @@ private static float[] lineConvolve(float[] arr, float[] kernel, boolean flip) { int c = 0; for (int k = -kw; k <= kw; k++) { int q = i - k; - if (q < arr.length && 0 <= q) { - y[i] += arr[q] * kernel[c]; - } else { - y[i] += arr[arr.length - 1] * kernel[c]; - } + if (q < arr.length && 0 <= q) y[i] += arr[q] * kernel[c]; + else y[i] += arr[arr.length - 1] * kernel[c]; c++; } } + return y; } @@ -502,9 +584,7 @@ public static void contrastAdjust(FloatProcessor fpaux, double dr, final double for (int i = 0; i < pixels.length; i++) { final int x = i % width; final int y = i / width; - if (rect.contains(x, y)) { - pixels[i] = (float) (pixels[i] * dr + d1); - } + if (rect.contains(x, y)) pixels[i] = (float) (pixels[i] * dr + d1); } } @@ -514,19 +594,18 @@ public static float[] findMinAndMax(FloatProcessor fp) { Rectangle rect = fp.getRoi(); float min = pixels[0]; float max = min; + for (int i = 0; i < pixels.length; i++) { final int x = i % width; final int y = i / width; if (rect.contains(x, y)) { float value = pixels[i]; if (!Float.isInfinite(value)) { - if (value < min) - min = value; - if (value > max) - max = value; + if (value < min) min = value; + if (value > max) max = value; } } } return new float[]{min, max}; } -} \ No newline at end of file +}