kernels.cl.h

Go to the documentation of this file.

/*=========================================================================
This file is part of CustusX, an Image Guided Therapy Application.

Copyright (c) 2008-2014, SINTEF Department of Medical Technology
All rights reserved.

Redistribution and use in source and binary forms, with or without 
modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, 
   this list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above copyright notice, 
   this list of conditions and the following disclaimer in the documentation 
   and/or other materials provided with the distribution.

3. Neither the name of the copyright holder nor the names of its contributors 
   may be used to endorse or promote products derived from this software 
   without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE 
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR 
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER 
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, 
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=========================================================================*/

#ifndef KERNELS_CLH_
#define KERNELS_CLH_

/*******************/
/* Begin compile time definitions */
/*******************/
//MAX_PLANES - number of planes to include in closest planes
//N_PLANES - number of planes/images (the z dimension)
//METHOD - the reconstruction method
//PLANE_METHOD - the method used to find the closes planes
//MAX_MULTISTART_STARTS - how many guesses the algorithm should produce as starting points for finding the closes planes
//NEWNESS_FACTOR - Newness weight, should newer (pixels in newer planes) be preferred
//BRIGHTNESS_FACTOR - Brightness weight, when selecting pixrls, should similarity in intensity count

/*******************/
/* End compile time definitions */
/*******************/

/*******************/
/* Begin constants */
/*******************/

#define CUBE_SIZE 4

#define N_BLOCKS 10

// Reconstruction methods
#define METHOD_VNN 0
#define METHOD_VNN2 1
#define METHOD_DW 2
#define METHOD_ANISOTROPIC 3

// Plane searching methods
#define PLANE_HEURISTIC 0
#define PLANE_CLOSEST 1

// Anisotropic method specific constants
#define ANISOTROPIC_SIZE 3
#define ANISOTROPIC_GAUSSIAN_SIZE 3

#define ANISOTROPIC_METHOD_INCLUDE_ALL 0
#define ANISOTROPIC_METHOD_BILINEAR_ON_PLANE 1
#define ANISOTROPIC_METHOD_GAUSSIAN_ON_PLANE 2

#define ANISOTROPIC_WEIGHT_METHOD_DISTANCE 0
#define ANISOTROPIC_WEIGHT_METHOD_BRIGHTNESS 1
#define ANISOTROPIC_WEIGHT_METHOD_LATENESS 2
#define ANISOTROPIC_WEIGHT_METHOD_BOTH 3

#define ANISOTROPIC_METHOD ANISOTROPIC_METHOD_BILINEAR_ON_PLANE

// Local search related parameters

//#define LOCAL_SEARCH_DISTANCE (N_PLANES/(MAX_MULTISTART_STARTS)) 
//#define LOCAL_SEARCH_DISTANCE 10
/*****************/
/* End constants */
/*****************/

/*****************/
/* Begin structs */
/*****************/

typedef struct _close_plane
{
        float dist;
        short plane_id;
        unsigned char intensity;
        unsigned char padding; // Align with 4
} close_plane_t;

typedef struct _output_volume_type
{
        int3 size;
        float3 spacing;
        __global unsigned char* volume;
} output_volume_type;

/***************/
/* End structs */
/***************/

/****************/
/* Begin macros */
/****************/

//#define DEBUG
#define CHECK_PLANE_INDICES

#ifdef DEBUG
        #define DEBUG_PRINTF(...) if((get_global_id(0) % 5000) == 0) printf(##__VA_ARGS__)
        //#define DEBUG_PRINTF(...) printf(##__VA_ARGS__)
        //#define BOUNDS_CHECK(x, min, ma x) if(x < min || x >= max) printf("Line %d: %s out of range: %d min: %d max: %d\n", __LINE__, #x, x, min, max)
        #define BOUNDS_CHECK(x, min, max)
#else
//      #define DEBUG_PRINTF(...)
        #define BOUNDS_CHECK(x, min, max)
#endif

//#define PLANE_DIST(voxel, matrix) (dot(matrix.s26AE,voxel) - dot(matrix.s26AE, matrix.s37BF))

#define EUCLID_DIST(a, b, c) sqrt((a)*(a) + (b)*(b) + (c)*(c))

#define PROJECTONTOPLANE(voxel, matrix, dist) (voxel - dist*(matrix.s26AE))

#define PROJECTONTOPLANEEQ(voxel, eq, dist) (voxel - dist*(eq))

#define ISINSIDE(x, size) ((x) >= 0 && (x) < (size))
#define ISNOTMASKED(x, y, mask, xsize) ((mask)[(x) + (y)*(xsize)] > 0)
//#define ISNOTMASKED(x, y, mask, xsize) true

#define VOXEL(v,x,y,z) v[x + y*volume_xsize + z*volume_ysize*volume_xsize]

#define WEIGHT_INV(x) (1.0f/fabs(x))
#define WEIGHT_INV2(x) (1.0f/fabs(x*x))
#define WEIGHT_INV4(x) (1.0f/fabs(x*x*x*x))
#define WEIGHT_SUB(x) (1.0f - fabs(x))

#define WEIGHT_TERNARY(val, mean, factor)                            \
    ((val) >= (mean) ? (factor) : 0.0f)

#define ANISOTROPIC_GAUSS_WEIGHT(px, var, mean, mean_id, sigma) WEIGHT_GAUSS(px.dist, sigma)

#ifndef ANISOTROPIC_WEIGHT_METHOD
        #define ANISOTROPIC_WEIGHT_METHOD ANISOTROPIC_WEIGHT_METHOD_BOTH
#endif

#ifndef BRIGHTNESS_FACTOR
        #define BRIGHTNESS_FACTOR 5.0f
#endif

#ifndef NEWNESS_FACTOR
        #define NEWNESS_FACTOR 5.0f
#endif

#define ANISOTROPIC_WEIGHT_BRIGHTNESS(px, var, mean, mean_id, sigma)    \
    ((WEIGHT_GAUSS(px.dist, sigma)) + (WEIGHT_TERNARY(px.intensity, mean, BRIGHTNESS_FACTOR)))

#define ANISOTROPIC_WEIGHT_LATENESS(px, var, mean, mean_id, sigma)      \
    ((WEIGHT_GAUSS(px.dist, sigma)) + (WEIGHT_TERNARY(px.plane_id, mean_id, NEWNESS_FACTOR)))

#define ANISOTROPIC_WEIGHT_BOTH(px, var, mean, mean_id, sigma) \
    ((WEIGHT_GAUSS(px.dist, sigma)) \
     + (WEIGHT_TERNARY(px.plane_id, mean_id, NEWNESS_FACTOR)) \
     + (WEIGHT_TERNARY(px.intensity, mean, BRIGHTNESS_FACTOR)))

#if ANISOTROPIC_WEIGHT_METHOD == ANISOTROPIC_WEIGHT_METHOD_DISTANCE
        #define ANISOTROPIC_WEIGHT(px, var, mean, mean_id, sigma) ANISOTROPIC_GAUSS_WEIGHT(px, var, mean, mean_id, sigma)
#elif ANISOTROPIC_WEIGHT_METHOD == ANISOTROPIC_WEIGHT_METHOD_BRIGHTNESS
        #define ANISOTROPIC_WEIGHT(px, var, mean, mean_id, sigma) ANISOTROPIC_WEIGHT_BRIGHTNESS(px, var, mean, mean_id, sigma)
#elif ANISOTROPIC_WEIGHT_METHOD == ANISOTROPIC_WEIGHT_METHOD_LATENESS
        #define ANISOTROPIC_WEIGHT(px, var, mean, mean_id, sigma) ANISOTROPIC_WEIGHT_LATENESS(px, var, mean, mean_id, sigma)
#elif ANISOTROPIC_WEIGHT_METHOD == ANISOTROPIC_WEIGHT_METHOD_BOTH
        #define ANISOTROPIC_WEIGHT(px, var, mean, mean_id, sigma) ANISOTROPIC_WEIGHT_BOTH(px, var, mean, mean_id, sigma)
#endif

// Gaussian weight function
#define WEIGHT_GAUSS_SIGMA (0.05f)

#define WEIGHT_GAUSS_SQRT_2PI 2.506628275f

#define WEIGHT_GAUSS_NONEXP_PART(sigma) (1.0f/((sigma)*WEIGHT_GAUSS_SQRT_2PI))
#define WEIGHT_GAUSS_EXP_PART(dist, sigma) exp(-((dist)*(dist))/(2*(sigma)*(sigma)))

#define WEIGHT_GAUSS(x, sigma) (WEIGHT_GAUSS_NONEXP_PART(sigma)*WEIGHT_GAUSS_EXP_PART(x, sigma))

#define DW_WEIGHT(x) WEIGHT_INV(x)
#define VNN2_WEIGHT(x) WEIGHT_INV(x)

#define CLOSE_PLANE_IDX(p, i) p[get_local_id(0)*(MAX_PLANES+1)+(i)]

#define CREATE_OUTPUT_VOLUME_TYPE(name, in_size, in_spacing, in_volume) \
        output_volume_type name; \
        name.size = in_size; \
        name.spacing = in_spacing; \
        name.volume = in_volume;

/**************/
/* End macros */
/**************/

/********************/
/* Begin prototypes */
/********************/

// Declare all the functions, as Apple seems to need that
//---------------------DEBUGGING-FUNCTIONALITY---------------------
#ifdef DEBUG
        void printMatrix(float16 matrix);
#endif /* DEBUG */

//---------------------DEBUGGING-FUNCTIONALITY---------------------

int isValidPixel(int2 point, const __global unsigned char* mask, int2 in_size);

int findHighestIdx(__local close_plane_t *planes, int n);

int2 findClosestPlanes_heuristic(__local close_plane_t *close_planes,
                __local float4* const plane_eqs,
                __global float16* const plane_matrices,
                const float4 voxel,
                const float radius,
                int guess,
                bool doTermDistance,
                __global const unsigned char* mask,
                int2 in_size,
                float2 in_spacing);

int2 findClosestPlanes_multistart(__local close_plane_t *close_planes,
                __local float4* const plane_eqs,
                __global float16* const plane_matrices,
                const float4 voxel,
                const float radius,
                int *multistart_guesses,
                int n_multistart_guesses,
                bool doTermDistance,
                __global const unsigned char* mask,
                int2 in_size,
                float2 in_spacing);

#if PLANE_METHOD == PLANE_EXACT
        #define FIND_CLOSE_PLANES(a, b, c, d, e, f, g, h, i, j) findClosestPlanes_multistart(a, b, c, d, e, f, g, 1, h, i, j)
#elif PLANE_METHOD == PLANE_CLOSEST
        #ifdef MAX_MULTISTART_STARTS
        #undef MAX_MULTISTART_STARTS
        #define MAX_MULTISTART_STARTS 1
        #endif

        #define FIND_CLOSE_PLANES(a, b, c, d, e, f, g, h, i, j) findClosestPlanes_multistart(a, b, c, d, e, f, g, 0, h, i, j)
#endif

__global const unsigned char* getImageData(int plane_id,
                __global const unsigned char* bscans_blocks[],
                int2 in_size);

float4 transform(float16 matrix, float4 voxel);

float4 transform_inv(float16 matrix, float4 voxel);

float2 transform_inv_xy(float16 matrix, float4 voxel);

int2 round_int(float2 value);

int2 toImgCoord_int(float4 voxel, float16 plane_matrix, float2 in_spacing);

float2 toImgCoord_float(float4 voxel, float16 plane_matrix, float2 in_spacing);

float bilinearInterpolation(float x,
                float y,
                const __global unsigned char* image,
                int in_xsize);

#if METHOD == METHOD_VNN
        #define PERFORM_INTERPOLATION(a, b, c, d, e, f, g, h, i)          \
                        performInterpolation_vnn(a, b, c, d, e, f, g, h, i)
#elif METHOD == METHOD_VNN2
        #define PERFORM_INTERPOLATION(a, b, c, d, e, f, g, h, i)          \
                        performInterpolation_vnn2(a, b, c, d, e, f, g, h, i)
#elif METHOD == METHOD_DW
        #define PERFORM_INTERPOLATION(a, b, c, d, e, f, g, h, i)          \
                        performInterpolation_dw(a, b, c, d, e, f, g, h, i)
#elif METHOD == METHOD_ANISOTROPIC
        #define PERFORM_INTERPOLATION(a, b, c, d, e, f, g, h, i)          \
                        performInterpolation_anisotropic(a, b, c, d, e, f, g, h, i)
#endif

unsigned char performInterpolation_vnn(__local close_plane_t *close_planes,
                int n_close_planes,
                __global const float16 *plane_matrices,
                __local const float4 *plane_eqs,
                __global const unsigned char* bscans_blocks[],
                int2 in_size,
                float2 in_spacing,
                __global const unsigned char* mask,
                float4 voxel);

unsigned char performInterpolation_vnn2(__local close_plane_t *close_planes,
                int n_close_planes,
                __global const float16 *plane_matrices,
                __local const float4 *plane_eqs,
                __global const unsigned char* bscans_blocks[],
                int2 in_size,
                float2 in_spacing,
                __global const unsigned char* mask,
                float4 voxel);

unsigned char performInterpolation_dw(__local close_plane_t *close_planes,
                int n_close_planes,
                __global const float16 *plane_matrices,
                __local const float4 *plane_eqs,
                __global const unsigned char* bscans_blocks[],
                int2 in_size,
                float2 in_spacing,
                __global const unsigned char* mask,
                float4 voxel);

unsigned char performInterpolation_anisotropic(__local close_plane_t *close_planes,
                int n_close_planes,
                __global const float16 *plane_matrices,
                __local const float4 *plane_eqs,
                __global const unsigned char* bscans_blocks[],
                int2 in_size,
                float2 in_spacing,
                __global const unsigned char* mask,
                float4 voxel);

unsigned char anisotropicFilter(__local const close_plane_t *pixels,
                int n_planes);

void prepare_plane_eqs(__global float16 *plane_matrices,
                __local float4 *plane_eqs);

int findLocalMinimas(int *guesses,
                __local const float4 *plane_eqs,
                float radius,
                float4 voxel,
                float3 out_spacing,
                __global const float16 *plane_matrices,
                __global const unsigned char *mask);

__kernel void voxel_methods(int volume_xsize,
                int volume_ysize,
                int volume_zsize,
                float volume_xspacing,
                float volume_yspacing,
                float volume_zspacing,
                int in_xsize,
                int in_ysize,
                float in_xspacing,
                float in_yspacing,
                // TODO: Wouldn't it be kind of nice if the bscans was an image sampler object?
                __global unsigned char* in_bscans_b0,
                __global unsigned char* in_bscans_b1,
                __global unsigned char* in_bscans_b2,
                __global unsigned char* in_bscans_b3,
                __global unsigned char* in_bscans_b4,
                __global unsigned char* in_bscans_b5,
                __global unsigned char* in_bscans_b6,
                __global unsigned char* in_bscans_b7,
                __global unsigned char* in_bscans_b8,
                __global unsigned char* in_bscans_b9,
                __global unsigned char* out_volume,
                __global float16 *plane_matrices,
                __global unsigned char* mask,
                __local float4 *plane_eqs,
                __local close_plane_t *planes,
                float radius);

/******************/
/* End prototypes */
/******************/

#endif /* KERNELS_CLH_ */