TomographicImaging/CCPi-Regularisation-Toolkit

Spatially variant regularisation field

dkazanc opened this issue · 10 comments

dkazanc commented
  • One requires a capability of passing not only a scalar but a vector of regularisation values.
  • When the N-dimensional field of regularisation parameters is evaluated, it can be a very powerful tool to incorporate better models into an algorithm making it robust against data inconsistencies.

Current implementation accepts only scalars for regularisation parameter so the core must be modified.

gfardell commented

I've had a quick look at this function.

A simple way without changing the code too much, and avoiding repetition would be to pass a pointer to lambda and a flag on whether it's an array or not. Then keeping your code mostly the same you just need to step through the array if it's an array, and don't step if it's not:

So here lambda_is_arr equals 1 or 0.

and the value used in your final calculation is then:

lambda_val = *(lambda + index* lambda_is_arr);

It's not a very safe solution - you'll have to make check the size of lambda beforehand.

So an example would be something like this:

float TV_kernel(float *D1, float *D2, float *D3, float *B, float *A, float * lambda, int lambda_is_arr, float tau, long dimX, long dimY, long dimZ)
{
	float  lambda_val;

	float dv1, dv2, dv3;
	long index, i, j, k, i1, i2, k1, j1, j2, k2;

	if (dimZ > 1)
	{
#pragma omp parallel for shared (D1, D2, D3, B, dimX, dimY, dimZ) private(index, i, j, k, i1, j1, k1, i2, j2, k2, dv1,dv2,dv3,lambda_val)
		for (j = 0; j < dimY; j++)
		{
			for (i = 0; i < dimX; i++)
			{
				for (k = 0; k < dimZ; k++)
				{

					index = (dimX*dimY)*k + j * dimX + i;

					lambda_val = *(lambda + index* lambda_is_arr);

					/* symmetric boundary conditions (Neuman) */
					i1 = i + 1; if (i1 >= dimX) i1 = i - 1;
					i2 = i - 1; if (i2 < 0) i2 = i + 1;
					j1 = j + 1; if (j1 >= dimY) j1 = j - 1;
					j2 = j - 1; if (j2 < 0) j2 = j + 1;
					k1 = k + 1; if (k1 >= dimZ) k1 = k - 1;
					k2 = k - 1; if (k2 < 0) k2 = k + 1;

					/*divergence components */
					dv1 = D1[index] - D1[(dimX*dimY)*k + j2 * dimX + i];
					dv2 = D2[index] - D2[(dimX*dimY)*k + j * dimX + i2];
					dv3 = D3[index] - D3[(dimX*dimY)*k2 + j * dimX + i];

					B[index] += tau * (lambda_val*(dv1 + dv2 + dv3) - (B[index] - A[index]));
				}
			}
		}
	}
	else
	{
#pragma omp parallel for shared (D1, D2, B, dimX, dimY) private(index, i, j, i1, j1, i2, j2,dv1,dv2,lambda_val)
		for (j = 0; j < dimY; j++)
		{
			for (i = 0; i < dimX; i++)
			{
				index = j * dimX + i;
				lambda_val = *(lambda + index * lambda_is_arr);

				/* symmetric boundary conditions (Neuman) */
				i1 = i + 1; if (i1 >= dimX) i1 = i - 1;
				i2 = i - 1; if (i2 < 0) i2 = i + 1;
				j1 = j + 1; if (j1 >= dimY) j1 = j - 1;
				j2 = j - 1; if (j2 < 0) j2 = j + 1;

				/* divergence components  */
				dv1 = D1[index] - D1[j2*dimX + i];
				dv2 = D2[index] - D2[j*dimX + i2];

				B[index] += tau * (lambda_val*(dv1 + dv2) - (B[index] - A[index]));
			}
		}
	}
	return *B;
}
dkazanc commented

Thanks a lot @gfardell, I find this as a quite neat solution actually! The problem I'm currently having is how to process inputs in Cython.
E.g. this function/wrapper

def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                     np.ndarray[np.float32_t, ndim=2, mode="c"] regularisation_parameter,
                     int iterationsNumb,
                     float marching_step_parameter,
                     float tolerance_param):

So the regularisation_parameter is defined and fixed as a 2D array, i.e. I cannot pass a scalar to it. Is there any other way around it? cheers

dkazanc commented

Actually I think I know what I should do. I need to check what regularisation_parameter is before calling the function.

def TV_ROF_CPU(inputData, regularisation_parameter, iterationsNumb, marching_step_parameter,tolerance_param):
    if inputData.ndim == 2:
      #### check regularisation_parameter
        return TV_ROF_2D(inputData, regularisation_parameter, iterationsNumb, marching_step_parameter,tolerance_param)
dkazanc commented

probably not very neat but this is what I do:

def TV_ROF_CPU(inputData, regularisation_parameter, iterationsNumb, marching_step_parameter,tolerance_param):
    if inputData.ndim == 2:
        if regularisation_parameter.ndim == 2:
            return TV_ROF_2D_vecreg(inputData, regularisation_parameter, iterationsNumb, marching_step_parameter,tolerance_param)
        else:
            return TV_ROF_2D_scalreg(inputData, regularisation_parameter, iterationsNumb, marching_step_parameter,tolerance_param)

Then there are 2 separate functions dealing with scalar or vector regularisation_parameter

paskino commented

I suppose you are looking at this

I'd probably change it to something like this

def TV_ROF_CPU(inputData, regularisation_parameter, iterationsNumb, marching_step_parameter,tolerance_param):
    cdef float lambdareg;
    if inputData.ndim == 2:
        if isinstance (regularisation_parameter, np.ndarray):
            return TV_ROF_2D(inputData, &regularisation_parameter[0,0], 1,  iterationsNumb, marching_step_parameter,tolerance_param)
        else: # supposedly this would be a float
            lambdareg = regularisation_parameter
            return TV_ROF_2D(inputData, &lambdareg, 0,  iterationsNumb, marching_step_parameter,tolerance_param)
    elif inputData.ndim == 3:
        return TV_ROF_3D(inputData, regularisation_parameter, iterationsNumb, marching_step_parameter,tolerance_param)

def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                     float * regularisation_parameter, int is_array, 
                     int iterationsNumb,
                     float marching_step_parameter,
                     float tolerance_param):        
cdef long dims[2]
    dims[0] = inputData.shape[0]
    dims[1] = inputData.shape[1]

    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
            np.zeros([dims[0],dims[1]], dtype='float32')
    cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
            np.ones([2], dtype='float32')

    # Run ROF iterations for 2D data
    TV_ROF_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], regularisation_parameter, is_array, iterationsNumb, marching_step_parameter, tolerance_param, dims[1], dims[0], 1)

    return (outputData,infovec)

EDITED

dkazanc commented

Thanks. I've got 2 errors:

return TV_ROF_2D(inputData, &regularisation_parameter[0,0], 1,  iterationsNumb, marching_step_parameter,tolerance_param)

src/cpu_regularisers.pyx:45:40: Cannot take address of Python object```

and

return TV_ROF_2D(inputData, &lambdareg, 0, iterationsNumb, marching_step_parameter,tolerance_param)
src/cpu_regularisers.pyx:48:40: Cannot convert 'float *' to Python object

both errors are about the second parameter in TV_ROF_2D

paskino commented

I suppose that's right. We need to do that check in the python function and pass the C function what it expects. Try this

def TV_ROF_CPU(inputData, regularisation_parameter, iterationsNumb, marching_step_parameter,tolerance_param):

    if inputData.ndim == 2:
        return TV_ROF_2D(inputData, regularisation_parameter,  iterationsNumb, marching_step_parameter,tolerance_param)
    elif inputData.ndim == 3:
        return TV_ROF_3D(inputData, regularisation_parameter, iterationsNumb, marching_step_parameter,tolerance_param)

def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                     regularisation_parameter, 
                     int iterationsNumb,
                     float marching_step_parameter,
                     float tolerance_param):        
    cdef long dims[2]
    dims[0] = inputData.shape[0]
    dims[1] = inputData.shape[1]
    cdef float lambdareg;
    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
            np.zeros([dims[0],dims[1]], dtype='float32')
    cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
            np.ones([2], dtype='float32')
    # Run ROF iterations for 2D data
    if isinstance (regularisation_parameter, np.ndarray):
        TV_ROF_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], &regularisation_parameter[0,0],  1, iterationsNumb, marching_step_parameter, tolerance_param, dims[1], dims[0], 1)

    else: # supposedly this would be a float
        lambdareg = regularisation_parameter;
        TV_ROF_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], &lambdareg,  0, iterationsNumb, marching_step_parameter, tolerance_param, dims[1], dims[0], 1)


    return (outputData,infovec)
dkazanc commented

thanks, we still get the Cannot take address of Python object error for &regularisation_parameter[0,0]

paskino commented

OK, the reason is, I presume, that regularisation_parameter is a Python object, whatever that may be.
I think one solution is to copy the regularisation_parameter array to a cdef array.

I haven't tested this

def TV_ROF_CPU(inputData, regularisation_parameter, iterationsNumb, marching_step_parameter,tolerance_param):

    if inputData.ndim == 2:
        return TV_ROF_2D(inputData, regularisation_parameter,  iterationsNumb, marching_step_parameter,tolerance_param)
    elif inputData.ndim == 3:
        return TV_ROF_3D(inputData, regularisation_parameter, iterationsNumb, marching_step_parameter,tolerance_param)

def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                     regularisation_parameter, 
                     int iterationsNumb,
                     float marching_step_parameter,
                     float tolerance_param):        
    cdef long dims[2]
    dims[0] = inputData.shape[0]
    dims[1] = inputData.shape[1]
    cdef float lambdareg;
    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] reg;
    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
            np.zeros([dims[0],dims[1]], dtype='float32')
    cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
            np.ones([2], dtype='float32')
    # Run ROF iterations for 2D data
    if isinstance (regularisation_parameter, np.ndarray):
        reg = regularisation_parameter.copy()
        TV_ROF_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], &reg[0,0],  1, iterationsNumb, marching_step_parameter, tolerance_param, dims[1], dims[0], 1)

    else: # supposedly this would be a float
        lambdareg = regularisation_parameter;
        TV_ROF_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], &lambdareg,  0, iterationsNumb, marching_step_parameter, tolerance_param, dims[1], dims[0], 1)


    return (outputData,infovec)
dkazanc commented

OK, as if by magic - all works! Thanks @gfardell and @paskino ! I'm already removing artifacts in reconstruction with that. Will adapt all the modules tmr and do PR.
reg_var