changing direct conversion to CDP, handling sparse tasks (3c30a0bc) · Commits · Elphel / tile_processor_gpu

src/TileProcessor.cuh

+215 −7

Original line number	Diff line number	Diff line
		@@ -104,6 +104,8 @@ GPU run time =523.451927ms, (direct conversion: 24.080189999999998ms, imclt: 17.
		#define KERNELS_STEP (1 << KERNELS_LSTEP)
		#define TILESX (IMG_WIDTH / DTT_SIZE)
		#define TILESY (IMG_HEIGHT / DTT_SIZE)
		#define CONVERT_DIRECT_INDEXING_THREADS_LOG2 5
		#define CONVERT_DIRECT_INDEXING_THREADS (1 << CONVERT_DIRECT_INDEXING_THREADS_LOG2) // 32
		// Make TILESYA >= TILESX and a multiple of 4
		#define TILESYA ((TILESY +3) & (~3))
		// increase row length by 1 so vertical passes will use different ports
		@@ -1541,10 +1543,132 @@ __global__ void gen_texture_list(

		#endif //#ifdef USE_CDP

		//#define CONVERT_DIRECT_INDEXING_THREADS_LOG2 5
		//#define CONVERT_DIRECT_INDEXING_THREADS (1 << CONVERT_DIRECT_INDEXING_THREADS_LOG2) // 32
		//#define CONVERT_DIRECT_NUM_CHUNKS ((TILESY*TILESX+CONVERT_DIRECT_INDEXING_THREADS-1) >> CONVERT_DIRECT_INDEXING_THREADS_LOG2)
		//#define CONVERT_DIRECT_NUM_CHUNKS2 ((CONVERT_DIRECT_NUM_CHUNKS+CONVERT_DIRECT_INDEXING_THREADS-1) >> CONVERT_DIRECT_INDEXING_THREADS_LOG2)
		//__global__ int num_active_tiles;
		//__global__ int active_tiles [TILESY*TILESX]; // indices of tiles in gpu_tasks that have non-zero correlations and/or textures
		//__device__ int num_acive_per_chunk [CONVERT_DIRECT_NUM_CHUNKS+1];
		//__device__ int num_acive_per_chunk2[CONVERT_DIRECT_NUM_CHUNKS2+1];

		// not maintaining order of the tiles to be processed
		extern "C" __global__ void index_direct(
		struct tp_task * gpu_tasks,
		int num_tiles, // number of tiles in task
		int * active_tiles, // pointer to the calculated number of non-zero tiles
		int * num_active_tiles) // indices to gpu_tasks // should be initialized to zero
		{
		int num_tile = blockIdx.x * blockDim.x + threadIdx.x;
		if (num_tile >= num_tiles){
		return;
		}
		if (gpu_tasks[num_tile].task != 0) {
		active_tiles[atomicAdd(num_active_tiles, 1)] = num_tile;
		}
		}

		extern "C" __global__ void convert_direct( // called with a single block, CONVERT_DIRECT_INDEXING_THREADS threads
		// struct CltExtra ** gpu_kernel_offsets, // [NUM_CAMS], // changed for jcuda to avoid struct parameters
		float ** gpu_kernel_offsets, // [NUM_CAMS],
		float ** gpu_kernels, // [NUM_CAMS],
		float ** gpu_images, // [NUM_CAMS],
		struct tp_task * gpu_tasks,
		float ** gpu_clt, // [NUM_CAMS][TILESY][TILESX][NUM_COLORS][DTT_SIZE*DTT_SIZE]
		size_t dstride, // in floats (pixels)
		int num_tiles, // number of tiles in task
		int lpf_mask, // apply lpf to colors : bit 0 - red, bit 1 - blue, bit2 - green. Now - always 0 !
		int woi_width,
		int woi_height,
		int kernels_hor,
		int kernels_vert,
		int * gpu_active_tiles, // pointer to the calculated number of non-zero tiles
		int * pnum_active_tiles) // indices to gpu_tasks
		{
		dim3 threads0(CONVERT_DIRECT_INDEXING_THREADS, 1, 1);
		dim3 blocks0 ((num_tiles + CONVERT_DIRECT_INDEXING_THREADS -1) >> CONVERT_DIRECT_INDEXING_THREADS_LOG2,1, 1);
		if (threadIdx.x == 0) { // of CONVERT_DIRECT_INDEXING_THREADS
		*pnum_active_tiles = 0;
		index_direct<<<blocks0,threads0>>>(
		gpu_tasks, // struct tp_task * gpu_tasks,
		num_tiles, //int num_tiles, // number of tiles in task
		gpu_active_tiles, //int * active_tiles, // pointer to the calculated number of non-zero tiles
		pnum_active_tiles); //int * pnum_active_tiles) // indices to gpu_tasks // should be initialized to zero
		cudaDeviceSynchronize();
		// now call actual convert_correct_tiles
		dim3 threads_tp(THREADSX, TILES_PER_BLOCK, 1);
		dim3 grid_tp((*pnum_active_tiles + TILES_PER_BLOCK -1 )/TILES_PER_BLOCK, 1);
		convert_correct_tiles<<<grid_tp,threads_tp>>>(
		gpu_kernel_offsets, // float ** gpu_kernel_offsets, // [NUM_CAMS],
		gpu_kernels, // float ** gpu_kernels, // [NUM_CAMS],
		gpu_images, // float ** gpu_images, // [NUM_CAMS],
		gpu_tasks, // struct tp_task * gpu_tasks, // array of tasks
		gpu_active_tiles, // int * gpu_active_tiles, // indices in gpu_tasks to non-zero tiles
		*pnum_active_tiles, // int num_active_tiles, // number of tiles in task
		gpu_clt, // float ** gpu_clt, // [NUM_CAMS][TILESY][TILESX][NUM_COLORS][DTT_SIZE*DTT_SIZE]
		dstride, // size_t dstride, // in floats (pixels)
		lpf_mask, // int lpf_mask, // apply lpf to colors : bit 0 - red, bit 1 - blue, bit2 - green. Now - always 0 !
		woi_width, // int woi_width, // varaible to swict between EO and LWIR
		woi_height, // int woi_height, // varaible to swict between EO and LWIR
		kernels_hor, // int kernels_hor, // varaible to swict between EO and LWIR
		kernels_vert); // int kernels_vert); // varaible to swict between EO and LWIR
		}
		}
		#if 0 // trying to keep the same order
		extern "C" __global__ void index_direct_init(
		int num_chunks, // number of tiles in task
		struct convert_direct_tmp* tmp)
		{
		int chunk_index = blockIdx.x * blockDim.x + threadIdx.x;
		if (chunk_index <= num_chunks){
		tmp->num_acive_per_chunk[chunk_index] = 0;
		}
		}
		extern "C" __global__ void index_direct(
		struct tp_task * gpu_tasks,
		int num_tiles, // number of tiles in task
		struct convert_direct_tmp* tmp)
		{
		__shared__ int num_active;
		if (threadIdx.x == 0) {
		num_active = 0;
		}
		__syncthreads();
		int num_tile = blockIdx.x * blockDim.x + threadIdx.x;
		if (num_tile < num_tiles) {
		if (gpu_tasks[num_tile].task){
		atomicAdd(&num_active, 1);
		}
		}
		__syncthreads();
		tmp -> num_acive_per_chunk[(num_tile >> CONVERT_DIRECT_INDEXING_THREADS_LOG2) + 1] = num_active; // skip [0]
		}

		extern "C"
		__global__ void convert_correct_tiles(
		extern "C" __global__ void build_index_direct(
		struct tp_task * gpu_tasks,
		int num_tiles, // number of tiles in task
		struct convert_direct_tmp* tmp)
		{
		__shared__ int num_active;
		if (threadIdx.x == 0) {
		num_active = 0;
		}
		__syncthreads();
		int num_tile = blockIdx.x * blockDim.x + threadIdx.x;
		if (num_tile < num_tiles) {
		if (gpu_tasks[num_tile].task){
		atomicAdd(&num_active, 1);
		}
		}
		__syncthreads();
		tmp -> num_acive_per_chunk[(num_tile >> CONVERT_DIRECT_INDEXING_THREADS_LOG2) + 1] = num_active; // skip [0]
		}


		/**
		* Top level to call other kernel with CDP
		*/
		extern "C" __global__ void convert_direct( // called with a single block, CONVERT_DIRECT_INDEXING_THREADS threads
		// struct CltExtra ** gpu_kernel_offsets, // [NUM_CAMS], // changed for jcuda to avoid struct parameters
		float ** gpu_kernel_offsets, // [NUM_CAMS],
		float ** gpu_kernels, // [NUM_CAMS],
		@@ -1553,6 +1677,85 @@ __global__ void convert_correct_tiles(
		float ** gpu_clt, // [NUM_CAMS][TILESY][TILESX][NUM_COLORS][DTT_SIZE*DTT_SIZE]
		size_t dstride, // in floats (pixels)
		int num_tiles, // number of tiles in task
		int lpf_mask, // apply lpf to colors : bit 0 - red, bit 1 - blue, bit2 - green. Now - always 0 !
		int woi_width,
		int woi_height,
		int kernels_hor,
		int kernels_vert,
		int * num_active_tiles,
		struct convert_direct_tmp* tmp) // temporary storage - avoiding static data for future overlap of kernel execution
		{
		int num_chunks = (num_tiles + CONVERT_DIRECT_INDEXING_THREADS -1) >> CONVERT_DIRECT_INDEXING_THREADS_LOG2;
		int num_chunk_blocks = (num_chunks + CONVERT_DIRECT_INDEXING_THREADS -1) >> CONVERT_DIRECT_INDEXING_THREADS_LOG2;
		dim3 threads0(CONVERT_DIRECT_INDEXING_THREADS, 1, 1);
		dim3 blocks0 (num_chunk_blocks,1, 1);
		__shared__ int superchunks[CONVERT_DIRECT_INDEXING_THREADS + 1];
		if (threadIdx.x == 0) { // of CONVERT_DIRECT_INDEXING_THREADS

		index_direct_init<<<blocks0,threads0>>>(num_chunks, tmp); // zero num_acive_per_chunk[]

		cudaDeviceSynchronize(); // not needed yet, just for testing
		index_direct<<<blocks0,threads0>>>(
		gpu_tasks, // struct tp_task * gpu_tasks,
		num_tiles, // int num_tiles)
		tmp);
		cudaDeviceSynchronize(); // not needed yet, just for testing
		// single-threaded - make cumulative
		// tmp-> num_acive_per_chunk2[0] = 0;
		superchunks[0] = 0;
		}
		__syncthreads();

		// calculate cumulative in 3 steps
		//1. num_acive_per_chunk2 with each element being a sum of CONVERT_DIRECT_INDEXING_THREADS (32) elements of num_acive_per_chunk
		int num_passes = (num_chunk_blocks + CONVERT_DIRECT_INDEXING_THREADS - 1) >> CONVERT_DIRECT_INDEXING_THREADS_LOG2;
		for (int pass = 0; pass < num_passes; pass++){
		int num_cluster2 = (pass << CONVERT_DIRECT_INDEXING_THREADS_LOG2) + threadIdx.x + 1; // skip 0
		if (num_cluster2 <= num_chunk_blocks){
		superchunks[threadIdx.x+1] = superchunks[0];
		int indx = num_cluster2 << CONVERT_DIRECT_INDEXING_THREADS_LOG2 + 1;
		for (int i = 0; i < CONVERT_DIRECT_INDEXING_THREADS; i++){
		if (indx <= num_chunks) {
		superchunks[threadIdx.x+1] += tmp -> num_acive_per_chunk[indx++];
		}
		}
		}
		__syncthreads();
		// make superchunks cumulative (single-threaded
		if (threadIdx.x == 0) { // of CONVERT_DIRECT_INDEXING_THREADS
		for (int i = 0; i < CONVERT_DIRECT_INDEXING_THREADS; i++){
		superchunks[i + 1] += superchunks[i];
		}
		}
		__syncthreads();
		// now update tmp -> num_acive_per_chunk[] by adding them together and adding the initial value

		if (num_cluster2 <= num_chunk_blocks){
		int indx = num_cluster2 << CONVERT_DIRECT_INDEXING_THREADS_LOG2 + 1;
		tmp -> num_acive_per_chunk[indx] += superchunks[threadIdx.x];
		for (int i = 0; i < CONVERT_DIRECT_INDEXING_THREADS; i++){
		int prev = tmp -> num_acive_per_chunk[indx++];
		if (indx <= num_chunks) {
		tmp -> num_acive_per_chunk[indx] += prev;
		}
		}
		}
		__syncthreads();
		}
		}
		#endif


		extern "C" __global__ void convert_correct_tiles(
		float ** gpu_kernel_offsets, // [NUM_CAMS],
		float ** gpu_kernels, // [NUM_CAMS],
		float ** gpu_images, // [NUM_CAMS],
		struct tp_task * gpu_tasks,
		int * gpu_active_tiles, // indices in gpu_tasks to non-zero tiles
		int num_active_tiles, // number of tiles in task
		float ** gpu_clt, // [NUM_CAMS][TILESY][TILESX][NUM_COLORS][DTT_SIZE*DTT_SIZE]
		size_t dstride, // in floats (pixels)
		// int num_tiles, // number of tiles in task
		int lpf_mask, // apply lpf to colors : bit 0 - red, bit 1 - blue, bit2 - green. Now - always 0 !
		int woi_width,
		int woi_height,
		@@ -1561,8 +1764,14 @@ __global__ void convert_correct_tiles(
		{
		dim3 t = threadIdx;
		int tile_in_block = threadIdx.y;
		int task_num = blockIdx.x * TILES_PER_BLOCK + tile_in_block;
		if (task_num >= num_tiles) return; // nothing to do
		// int task_num = blockIdx.x * TILES_PER_BLOCK + tile_in_block;
		// if (task_num >= num_tiles) return; // nothing to do
		int task_indx = blockIdx.x * TILES_PER_BLOCK + tile_in_block;
		if (task_indx >= num_active_tiles){
		return; // nothing to do
		}
		int task_num = gpu_active_tiles[task_indx];

		struct tp_task * gpu_task = &gpu_tasks[task_num];
		if (!gpu_task->task) return; // NOP tile
		__shared__ struct tp_task tt [TILES_PER_BLOCK];
		@@ -1626,8 +1835,7 @@ __global__ void convert_correct_tiles(
		woi_height, // int woi_height,
		kernels_hor, // int kernels_hor,
		kernels_vert); //int kernels_vert)

		__syncthreads();// __syncwarp();
		__syncthreads();
		}
		}
		}
		@@ -1989,7 +2197,7 @@ __global__ void textures_accumulate(
		int tileX = tile_num - tileY * TILESX;
		int tile_x0 = (tileX - (woi + 0)) DTT_SIZE; // - (DTT_SIZE/2); // may be negative == -4
		int tile_y0 = (tileY - (woi + 1)) DTT_SIZE; // - (DTT_SIZE/2); // may be negative == -4
		int height = *(woi + 3) << DTT_SIZE_LOG2;
		/// int height = *(woi + 3) << DTT_SIZE_LOG2;

		#ifdef DEBUG12
		if ((tile_num == DBG_TILE) && (threadIdx.x == 0) && (threadIdx.y == 0)){

src/TileProcessor.h

+36 −13

Original line number	Diff line number	Diff line
		@@ -41,9 +41,14 @@
		#include "tp_defines.h"
		#endif

		extern "C" __global__ void index_direct(
		struct tp_task * gpu_tasks,
		int num_tiles, // number of tiles in task
		int * active_tiles, // pointer to the calculated number of non-zero tiles
		int * num_active_tiles); // indices to gpu_tasks // should be initialized to zero

		extern "C"
		__global__ void convert_correct_tiles(
		extern "C" __global__ void convert_direct( // called with a single block, CONVERT_DIRECT_INDEXING_THREADS threads
		// struct CltExtra ** gpu_kernel_offsets, // [NUM_CAMS], // changed for jcuda to avoid struct parameters
		float ** gpu_kernel_offsets, // [NUM_CAMS],
		float ** gpu_kernels, // [NUM_CAMS],
		float ** gpu_images, // [NUM_CAMS],
		@@ -51,6 +56,24 @@ __global__ void convert_correct_tiles(
		float ** gpu_clt, // [NUM_CAMS][TILESY][TILESX][NUM_COLORS][DTT_SIZE*DTT_SIZE]
		size_t dstride, // in floats (pixels)
		int num_tiles, // number of tiles in task
		int lpf_mask, // apply lpf to colors : bit 0 - red, bit 1 - blue, bit2 - green. Now - always 0 !
		int woi_width,
		int woi_height,
		int kernels_hor,
		int kernels_vert,
		int * gpu_active_tiles, // pointer to the calculated number of non-zero tiles
		int * pnum_active_tiles); // indices to gpu_tasks

		extern "C" __global__ void convert_correct_tiles(
		float ** gpu_kernel_offsets, // [NUM_CAMS],
		float ** gpu_kernels, // [NUM_CAMS],
		float ** gpu_images, // [NUM_CAMS],
		struct tp_task * gpu_tasks,
		int * gpu_active_tiles, // indices in gpu_tasks to non-zero tiles
		int num_active_tiles, // number of tiles in task
		float ** gpu_clt, // [NUM_CAMS][TILESY][TILESX][NUM_COLORS][DTT_SIZE*DTT_SIZE]
		size_t dstride, // in floats (pixels)
		// int num_tiles, // number of tiles in task
		int lpf_mask, // apply lpf to colors : bit 0 - red, bit 1 - blue, bit2 - green. Now - always 0 !
		int woi_width,
		int woi_height,

src/geometry_correction.cu

+65 −197

Original line number	Diff line number	Diff line
		@@ -62,6 +62,8 @@ __device__ void printExtrinsicCorrection(corr_vector * cv);
		inline __device__ float getRByRDist(float rDist,
		float rByRDist [RBYRDIST_LEN]); //shared memory



		__constant__ float ROTS_TEMPLATE[7][3][3][3] = {// ...{cos,sin,const}...
		{ // azimuth
		{{ 1, 0,0},{0, 0,0},{ 0,-1,0}},
		@@ -116,201 +118,6 @@ __constant__ int mm_seq [3][3][3]={
		{-1,-1,-1} // do nothing
		}};

		#if 0
		__device__ float rot_matrices [NUM_CAMS][3][3];
		//__device__ float rot_deriv_matrices [NUM_CAMS][4][3][3]; // /d_azimuth, /d_tilt, /d_roll, /d_zoom)


		// threads (3,3,4)
		extern "C" __global__ void calc_rot_matrices(
		struct corr_vector * gpu_correction_vector)
		{
		__shared__ float zoom [NUM_CAMS];
		__shared__ float sincos [NUM_CAMS][3][2]; // {az,tilt,roll, d_az, d_tilt, d_roll, d_az}{cos,sin}
		__shared__ float matrices[NUM_CAMS][4][3][3]; // [7] - extra

		float angle;
		int ncam = threadIdx.z;
		int nangle1 = threadIdx.x + threadIdx.y * blockDim.x; // * >> 1;
		int nangle = nangle1 >> 1;
		int is_sin = nangle1 & 1;

		#ifdef DEBUG20a
		if ((threadIdx.x == 0) && ( threadIdx.y == 0) && ( threadIdx.z == 0)){
		printf("\nget_tiles_offsets() threadIdx.x = %d, blockIdx.x= %d\n", (int)threadIdx.x, (int) blockIdx.x);
		printExtrinsicCorrection(gpu_correction_vector);
		}
		__syncthreads();// __syncwarp();
		#endif // DEBUG20


		if (nangle < 4){ // this part only for 1-st 3
		float* gangles =
		(nangle ==0)?gpu_correction_vector->azimuth:(
		(nangle ==1)?gpu_correction_vector->tilt:(
		(nangle ==2)?gpu_correction_vector->roll:
		gpu_correction_vector->zoom));
		if ((ncam < (NUM_CAMS -1)) \|\| (nangle == 2)){ // for rolls - all 4
		angle = *(gangles + ncam);

		} else {
		angle = 0.0f;

		#pragma unroll
		for (int n = 0; n < (NUM_CAMS-1); n++){
		angle -= *(gangles + n);
		}
		}
		if (!is_sin){
		angle += M_PI/2;
		}
		if (nangle < 3) {
		sincos[ncam][nangle][is_sin]=sinf(angle);
		} else if (is_sin){
		zoom[ncam] = angle;
		}
		}
		__syncthreads();


		#ifdef DEBUG20a
		if ((threadIdx.x == 0) && (threadIdx.y == 0) && (threadIdx.z == 0)){
		for (int n = 0; n < NUM_CAMS; n++){
		printf("\n Azimuth matrix for camera %d, sincos[0] = %f, sincos[1] = %f, zoom = %f\n", n, sincos[n][0][0], sincos[n][0][1], zoom[n]);
		printf(" Tilt matrix for camera %d, sincos[0] = %f, sincos[0] = %f\n", n, sincos[n][1][0], sincos[n][1][1]);
		printf(" Roll matrix for camera %d, sincos[0] = %f, sincos[2] = %f\n", n, sincos[n][2][0], sincos[n][2][1]);
		}
		}
		__syncthreads();// __syncwarp();
		#endif // DEBUG20


		if (nangle == 3) {
		sincos[ncam][2][is_sin] *= (1.0 + zoom[ncam]); // modify roll
		}
		__syncthreads();


		#ifdef DEBUG20a
		if ((threadIdx.x == 0) && (threadIdx.y == 0) && (threadIdx.z == 0)){
		for (int n = 0; n < NUM_CAMS; n++){
		printf("\na Azimuth matrix for camera %d, sincos[0] = %f, sincos[1] = %f, zoom = %f\n", n, sincos[n][0][0], sincos[n][0][1], zoom[n]);
		printf("a Tilt matrix for camera %d, sincos[0] = %f, sincos[0] = %f\n", n, sincos[n][1][0], sincos[n][1][1]);
		printf("a Roll matrix for camera %d, sincos[0] = %f, sincos[2] = %f\n", n, sincos[n][2][0], sincos[n][2][1]);
		}
		}
		__syncthreads();// __syncwarp();
		#endif // DEBUG20





		// now 3x3
		for (int axis = 0; axis < 3; axis++) {
		matrices[ncam][axis][threadIdx.y][threadIdx.x] =
		ROTS_TEMPLATE[axis][threadIdx.y][threadIdx.x][0] * sincos[ncam][axis][0]+ // cos
		ROTS_TEMPLATE[axis][threadIdx.y][threadIdx.x][1] * sincos[ncam][axis][1]+ // sin
		ROTS_TEMPLATE[axis][threadIdx.y][threadIdx.x][2]; // const
		}
		__syncthreads();


		#ifdef DEBUG20a
		if ((threadIdx.x == 0) && (threadIdx.y == 0) && (threadIdx.z == 0)){
		for (int n = 0; n < NUM_CAMS; n++){

		printf("\n1-Azimuth matrix for camera %d, sincos[0] = %f, sincos[1] = %f\n", n, sincos[n][0][0], sincos[n][0][1]);
		for (int i = 0; i < 3; i++){
		for (int j = 0; j < 3; j++){
		printf("%9.6f, ", matrices[n][0][i][j]);
		}
		printf("\n");
		}

		printf("1-Tilt matrix for camera %d, sincos[0] = %f, sincos[1] = %f\n", n, sincos[n][1][0], sincos[n][1][1]);
		for (int i = 0; i < 3; i++){
		for (int j = 0; j < 3; j++){
		printf("%9.6f, ", matrices[n][1][i][j]);
		}
		printf("\n");
		}

		printf("1-Roll/Zoom matrix for camera %d, sincos[0] = %f, sincos[1] = %f\n", n, sincos[n][2][0], sincos[n][2][1]);
		for (int i = 0; i < 3; i++){
		for (int j = 0; j < 3; j++){
		printf("%9.6f, ", matrices[n][2][i][j]);
		}
		printf("\n");
		}

		}
		}
		__syncthreads();// __syncwarp();
		#endif // DEBUG20





		// tilt * az ->
		// multiply matrices[ncam][1] * matrices[ncam][0] -> matrices[ncam][3]
		matrices[ncam][3][threadIdx.y][threadIdx.x] =
		matrices[ncam][1][threadIdx.y][0] * matrices[ncam][0][0][threadIdx.x]+
		matrices[ncam][1][threadIdx.y][1] * matrices[ncam][0][1][threadIdx.x]+
		matrices[ncam][1][threadIdx.y][2] * matrices[ncam][0][2][threadIdx.x];
		// multiply matrices[ncam][2] * matrices[ncam][3] -> rot_matrices[ncam]
		__syncthreads();
		rot_matrices[ncam][threadIdx.y][threadIdx.x] =
		matrices[ncam][2][threadIdx.y][0] * matrices[ncam][3][0][threadIdx.x]+
		matrices[ncam][2][threadIdx.y][1] * matrices[ncam][3][1][threadIdx.x]+
		matrices[ncam][2][threadIdx.y][2] * matrices[ncam][3][2][threadIdx.x];
		__syncthreads();


		#ifdef DEBUG20
		if ((threadIdx.x == 0) && (threadIdx.y == 0) && (threadIdx.z == 0)){
		for (int n = 0; n < NUM_CAMS; n++){

		printf("\n2 - Azimuth matrix for camera %d, sincos[0] = %f, sincos[1] = %f\n", n, sincos[n][0][0], sincos[n][0][1]);
		for (int i = 0; i < 3; i++){
		for (int j = 0; j < 3; j++){
		printf("%9.6f, ", matrices[n][0][i][j]);
		}
		printf("\n");
		}

		printf("2 - Tilt matrix for camera %d, sincos[0] = %f, sincos[1] = %f\n", n, sincos[n][1][0], sincos[n][1][1]);
		for (int i = 0; i < 3; i++){
		for (int j = 0; j < 3; j++){
		printf("%9.6f, ", matrices[n][1][i][j]);
		}
		printf("\n");
		}

		printf("2 - Roll/Zoom matrix for camera %d, sincos[0] = %f, sincos[1] = %f\n", n, sincos[n][2][0], sincos[n][2][1]);
		for (int i = 0; i < 3; i++){
		for (int j = 0; j < 3; j++){
		printf("%9.6f, ", matrices[n][2][i][j]);
		}
		printf("\n");
		}

		printf("2 - Rotation matrix for camera %d\n", n);
		for (int i = 0; i < 3; i++){
		for (int j = 0; j < 3; j++){
		printf("%9.6f, ", rot_matrices[n][i][j]);
		}
		printf("\n");
		}
		}
		}
		__syncthreads();// __syncwarp();
		#endif // DEBUG20


		}
		#endif
		__constant__ int offset_rots = 0; //0
		__constant__ int offset_derivs = 1; // 1..4 // should be next
		__constant__ int offset_matrices = 5; // 5..11
		@@ -890,8 +697,69 @@ extern "C" __global__ void get_tiles_offsets(


		}


		extern "C" __global__ void calcReverseDistortionTable(
		struct gc * geometry_correction,
		float * rByRDist)
		{
		//int num_threads = NUM_CAMS * blockDim.z * blockDim.y * blockDim.x; // 36
		int indx = ((blockIdx.x * blockDim.z + threadIdx.z) * blockDim.y + threadIdx.y) * blockDim.x + threadIdx.x;
		// double delta=1E-20; // 12; // 10; // -8; 215.983994 ms
		// double delta=1E-4; //rByRDist error = 0.000072
		double delta=1E-10; // 12; // 10; // -8; 0.730000 ms
		double minDerivative=0.01;
		int numIterations=1000;
		double drDistDr=1.0;
		double d=1.0
		-geometry_correction -> distortionA8
		-geometry_correction -> distortionA7
		-geometry_correction -> distortionA6
		-geometry_correction -> distortionA5
		-geometry_correction -> distortionA
		-geometry_correction -> distortionB
		-geometry_correction -> distortionC;
		double rPrev=0.0;
		int num_points = (RBYRDIST_LEN + CALC_REVERSE_TABLE_BLOCK_THREADS - 1) / CALC_REVERSE_TABLE_BLOCK_THREADS;
		for (int p = 0; p < num_points; p ++){
		int i = indx * num_points +p;
		if (i >= RBYRDIST_LEN){
		return;
		}
		if (i == 0){
		rByRDist[0]= (float) 1.0/d;
		break;
		}
		double rDist = RBYRDIST_STEP * i;
		double r = (p == 0) ? rDist : rPrev;
		for (int iteration=0;iteration<numIterations;iteration++){
		double k=(((((((
		geometry_correction -> distortionA8) * r +
		geometry_correction -> distortionA7) * r +
		geometry_correction -> distortionA6) * r +
		geometry_correction -> distortionA5) * r +
		geometry_correction -> distortionA) * r +
		geometry_correction -> distortionB) * r +
		geometry_correction -> distortionC) * r + d;
		drDistDr=(((((((
		8 * geometry_correction -> distortionA8) * r +
		7 * geometry_correction -> distortionA7) * r +
		6 * geometry_correction -> distortionA6) * r +
		5 * geometry_correction -> distortionA5) * r +
		4 * geometry_correction -> distortionA) * r +
		3 * geometry_correction -> distortionB) * r+
		2 * geometry_correction -> distortionC) * r+d;
		if (drDistDr<minDerivative) { // folds backwards !
		return; // too high distortion
		}
		double rD=r*k;
		if (fabs(rD-rDist)<delta){
		break;
		}
		r+=(rDist-rD)/drDistDr;
		}
		rPrev=r;
		rByRDist[i]= (float) r/rDist;
		}
		}

		/**
		* Calculate non-distorted radius from distorted using table approximation

src/geometry_correction.h

+6 −5

Original line number	Diff line number	Diff line
		@@ -148,14 +148,15 @@ extern "C" __global__ void get_tiles_offsets(
		float * gpu_rByRDist, // length should match RBYRDIST_LEN
		trot_deriv * gpu_rot_deriv);

		#if 0
		// uses 3 threadIdx.x, 3 - threadIdx.y, 4 - threadIdx.z
		extern "C" __global__ void calc_rot_matrices(
		struct corr_vector * gpu_correction_vector);
		#endif
		// uses NUM_CAMS blocks, (3,3,3) threads
		extern "C" __global__ void calc_rot_deriv(
		struct corr_vector * gpu_correction_vector,
		trot_deriv * gpu_rot_deriv);

		#define CALC_REVERSE_TABLE_BLOCK_THREADS (NUM_CAMS * 3 * 3 * 3) // fixed blockDim
		// Use same blocks/threads as with calc_rot_deriv() - NUM_CAMS blocks, (3,3,3) threads
		extern "C" __global__ void calcReverseDistortionTable(
		struct gc * geometry_correction,
		float * rByRDist);

src/test_tp.cu

+122 −51

File changed.

Preview size limit exceeded, changes collapsed.

Admin message