Is this algorithm applicable to 3D images of any width, height , and depth?

[Mechn]

I fix the code as follow:
`
// Input parameters
int width = 512;
int height = 512;
int depth = 256;
//int fboSize = 1024;
int nVertices = 100;

int phase1Band = 1;
const int phase2Band = 1;
int phase3Band = 1;

#define TOID(x, y, z, w, h) ((z) * (w) * (h) + (y) * (w) + (x))

// Generate input points
void generateRandomPoints(int width, int height, int depth, int nPoints)
{
int tx, ty, tz, id;

randinit(0);

for (int i = 0; i < width * height * depth; i++)
    inputVoronoi[i] = MARKER; 

for (int i = 0; i < nPoints; i++)
{
    do { 
        tx = int(random() * width);
        ty = int(random() * height);
        tz = int(random() * depth);
        id = TOID(tx, ty, tz, width, height); 
	} while (inputVoronoi[id] != MARKER); 

    inputVoronoi[id] = ENCODE(tx, ty, tz, 0, 0); 
    inputPoints[i] = ENCODE(tx, ty, tz, 0, 0);
}

}
// Initialization
void initialization()
{
pba3DInitialization(width, height, depth);

inputPoints     = (int *) malloc(nVertices * sizeof(int));
inputVoronoi    = (int *) malloc(width * height * depth * sizeof(int));
outputVoronoi   = (int *) malloc(width * height * depth * sizeof(int));

}
pha3DHost.cu
// Global Variables
int **pbaTextures;

size_t pbaMemSize;
int pbaCurrentBuffer;
//int pbaTexSize;
int pbaWidth;
int pbaHeight;
int pbaDepth;

// Initialize necessary memory for 3D Voronoi Diagram computation
// - textureSize: The size of the Discrete Voronoi Diagram (width = height)
void pba3DInitialization(int width, int height, int depth)
{
//pbaTexSize = fboSize;
pbaWidth = width;
pbaHeight = height;
pbaDepth = depth;

pbaTextures = (int **) std::malloc(2 * sizeof(int *)); 

pbaMemSize = width * height * depth * sizeof(int);

cudaMalloc((void **) &pbaTextures[0], pbaMemSize);
cudaMalloc((void **) &pbaTextures[1], pbaMemSize);

}

void pba3DColorZAxis(int m1)
{
dim3 block = dim3(BLOCKX, BLOCKY);
dim3 grid = dim3(pbaWidth / block.x, pbaHeight / block.y);

kernelFloodZ<<< grid, block >>>(pbaTextures[pbaCurrentBuffer], pbaTextures[1 - pbaCurrentBuffer], pbaWidth, pbaHeight, pbaDepth); 
pbaCurrentBuffer = 1 - pbaCurrentBuffer; 

}

void pba3DComputeProximatePointsYAxis(int m2)
{
dim3 block = dim3(BLOCKX, BLOCKY);
dim3 grid = dim3(pbaWidth / block.x, pbaDepth / block.y);

kernelMaurerAxis<<< grid, block >>>(pbaTextures[pbaCurrentBuffer], pbaTextures[1 - pbaCurrentBuffer], pbaWidth, pbaHeight, pbaDepth);

}

// Phase 3 of PBA. m3 must divides texture size
// This method color along the Y axis
void pba3DColorYAxis(int m3)
{
dim3 block = dim3(BLOCKSIZE, m3);
dim3 grid = dim3(pbaWidth / block.x, pbaHeight);

kernelColorAxis<<< grid, block >>>(pbaTextures[1 - pbaCurrentBuffer], pbaTextures[pbaCurrentBuffer], pbaWidth, pbaHeight, pbaDepth); 

}
pha3DKernel.h
// Flood along the Z axis
global void kernelFloodZ(int *input, int *output, int width, int height, int depth)
{
int tx = blockIdx.x * blockDim.x + threadIdx.x;
int ty = blockIdx.y * blockDim.y + threadIdx.y;
int tz = 0;

int plane = width * height; 
int id = TOID(tx, ty, tz, width, height); 
int pixel1, pixel2; 

pixel1 = ENCODE(0,0,0,1,0); 

// Sweep down
for (int i = 0; i < depth; i++, id += plane) {
    pixel2 = input[id];

    if (!NOTSITE(pixel2))
        pixel1 = pixel2;

    output[id] = pixel1;
}

int dist1, dist2, nz;

id -= plane + plane;

// Sweep up
for (int i = depth - 2; i >= 0; i--, id -= plane) {
    nz = GET_Z(pixel1);
    dist1 = abs(nz - (tz + i));

    pixel2 = output[id];
    nz = GET_Z(pixel2);
    dist2 = abs(nz - (tz + i));

    if (dist2 < dist1)
        pixel1 = pixel2;

    output[id] = pixel1;
}

}

global void kernelMaurerAxis(int *input, int *stack, int width, int height, int depth)
{
int tx = blockIdx.x * blockDim.x + threadIdx.x;
int tz = blockIdx.y * blockDim.y + threadIdx.y;
int ty = 0;

int id = TOID(tx, ty, tz, width, height);

int lasty = 0;
int x1, y1, z1, x2, y2, z2, nx, ny, nz;
int p = ENCODE(0,0,0,1,0), s1 = ENCODE(0,0,0,1,0), s2 = ENCODE(0,0,0,1,0);
int flag = 0;

for (ty = 0; ty < height; ++ty, id += width) {
    p = input[id];

    if (!NOTSITE(p)) {

        while (HASNEXT(s2)) {
            DECODE(s1, x1, y1, z1);
            DECODE(s2, x2, y2, z2);
            DECODE(p, nx, ny, nz);

            if (!dominate(x1, y2, z1, x2, lasty, z2, nx, ty, nz, tx, tz))
            	break;

            lasty = y2; s2 = s1; y2 = y1;

            if (HASNEXT(s2))
                s1 = stack[TOID(tx, y2, tz, width, height)];
        }

        DECODE(p, nx, ny, nz);
        s1 = s2;
        s2 = ENCODE(nx, lasty, nz, 0, flag);
        y2 = lasty;
        lasty = ty;

        stack[id] = s2;

        flag = 1;
    }
}

if (NOTSITE(p))
    stack[TOID(tx, ty - 1, tz, width, height)] = ENCODE(0, lasty, 0, 1, flag); 

}

global void kernelColorAxis(int *input, int *output, int width, int height, int depth)
{
shared int block[BLOCKSIZE][BLOCKSIZE];

int col = threadIdx.x;
int tid = threadIdx.y;
int tx = blockIdx.x * blockDim.x + col; 
int tz = blockIdx.y;

int x1, y1, z1, x2, y2, z2;
int last1 = ENCODE(0,0,0,1,0), last2 = ENCODE(0,0,0,1,0), lasty;
long long dx, dy, dz, best, dist;

lasty = height - 1;

last2 = input[TOID(tx, lasty, tz, width, height)]; 
DECODE(last2, x2, y2, z2);

if (NOTSITE(last2)) {
	lasty = y2;
	if(HASNEXT(last2)) {
		last2 = input[TOID(tx, lasty, tz, width, height)];
		DECODE(last2, x2, y2, z2);
	}
}

if (HASNEXT(last2)) {
	last1 = input[TOID(tx, y2, tz, width, height)];
	DECODE(last1, x1, y1, z1);
}

int y_start, y_end, n_step = width / blockDim.x;
for(int step = 0; step < n_step; ++step) {
	y_start = height - step * blockDim.x - 1;
	y_end = height - (step + 1) * blockDim.x;

    for (int ty = y_start - tid; ty >= y_end; ty -= blockDim.y) {
    	dx = x2 - tx; dy = lasty - ty; dz = z2 - tz;
		// Euclidean distance metric
		best = dx * dx + dy * dy + dz * dz;

		while (HASNEXT(last2)) {
			dx = x1 - tx; dy = y2 - ty; dz = z1 - tz;
			dist = dx * dx + dy * dy + dz * dz;

			if(dist > best) break;

			best = dist; lasty = y2; last2 = last1;
			DECODE(last2, x2, y2, z2);

			if (HASNEXT(last2)) {
				last1 = input[TOID(tx, y2, tz, width, height)];
				DECODE(last1, x1, y1, z1); 
			}
        }

        block[threadIdx.x][ty - y_end] = ENCODE(lasty, x2, z2, NOTSITE(last2), 0);
    }

    __syncthreads();

    if(!threadIdx.y) {
    	int id = TOID(y_end + threadIdx.x, blockIdx.x * blockDim.x, tz, width, height);
    	for(int i = 0; i < blockDim.x; i++, id+=width) {
    		output[id] = block[i][threadIdx.x];
    	}
    }

    __syncthreads();
}

}
`

Appreciate for answer! Thank you!!