Question about tvdenoising.cu

Hello. I have a question regarding the tvdenoising algorithm (CUDA/tvdenoising.cu file.).
Firstly, please forgive my lack of theoretical background knowledge on TV denoising.

My question is about the update_p function in the code below.

TIGRE/Common/CUDA/tvdenoising.cu

Lines 145 to 170 in ae39ee0

    
           __global__ 
        
                   void update_p(const float* u, float* pz, float* py, float* px, 
        
                   float tau, long depth, long rows, long cols, 
        
                   float dz, float dy, float dx) 
        
           { 
        
               long x = threadIdx.x + blockIdx.x * blockDim.x; 
        
               long y = threadIdx.y + blockIdx.y * blockDim.y; 
        
               long z = threadIdx.z + blockIdx.z * blockDim.z; 
        
               long idx = z * rows * cols + y * cols + x; 
        
               if ( x >= cols || y >= rows || z >= depth ) 
        
                   return; 
        
               float grad[3] = {0,0,0}, q[3]; 
        
               gradient(u, grad, z, y, x, depth, rows, cols, dz, dy, dx); 
        
               q[0] = pz[idx] + tau * grad[0]; 
        
               q[1] = py[idx] + tau * grad[1]; 
        
               q[2] = px[idx] + tau * grad[2]; 
        
               float norm = fmaxf(1.0f, sqrtf(q[0] * q[0] + q[1] * q[1] + q[2] * q[2])); 
        
               pz[idx] = q[0] / norm; 
        
               py[idx] = q[1] / norm; 
        
               px[idx] = q[2] / norm; 
        
           }

What does p represent here?

Initially, when looking at the update_u function, I thought p might be $\nabla u/ | \nabla u|$ from ROF PDE:

$$\nabla \cdot (\nabla u/ | \nabla u|) + \lambda(f-u) = 0$$

However, I'm having trouble understanding the logic of updating p in the update_p function.
Could you please explain this logic, or perhaps point me to a theoretical reference related to this?

Thank you for your time. I would really appreciate your assistance.

Its one of the variables for the primal-dual algorithm used to minimise the ROF.
The code is essentially an implementation of this paper: https://pubmed.ncbi.nlm.nih.gov/20352289/. You can read more in the paper or my PhD thesis.

I actually should add a reference to it, there is a dead link to the paper in the source code. Should update that.

Thank you very much. I'll check the paper!

	__global__
	void update_p(const float* u, float* pz, float* py, float* px,
	float tau, long depth, long rows, long cols,
	float dz, float dy, float dx)
	{
	long x = threadIdx.x + blockIdx.x * blockDim.x;
	long y = threadIdx.y + blockIdx.y * blockDim.y;
	long z = threadIdx.z + blockIdx.z * blockDim.z;
	long idx = z * rows * cols + y * cols + x;

	if ( x >= cols \|\| y >= rows \|\| z >= depth )
	return;

	float grad[3] = {0,0,0}, q[3];
	gradient(u, grad, z, y, x, depth, rows, cols, dz, dy, dx);

	q[0] = pz[idx] + tau * grad[0];
	q[1] = py[idx] + tau * grad[1];
	q[2] = px[idx] + tau * grad[2];

	float norm = fmaxf(1.0f, sqrtf(q[0] * q[0] + q[1] * q[1] + q[2] * q[2]));

	pz[idx] = q[0] / norm;
	py[idx] = q[1] / norm;
	px[idx] = q[2] / norm;
	}