Autonomous Space Robotics LabSpeeded Up SURF |
|
keypoint_interpolation.cu
Go to the documentation of this file.00001 /* 00002 Copyright (c) 2010, Paul Furgale and Chi Hay Tong 00003 All rights reserved. 00004 00005 Redistribution and use in source and binary forms, with or without 00006 modification, are permitted provided that the following conditions are 00007 met: 00008 00009 * Redistributions of source code must retain the above copyright notice, 00010 this list of conditions and the following disclaimer. 00011 * Redistributions in binary form must reproduce the above copyright 00012 notice, this list of conditions and the following disclaimer in the 00013 documentation and/or other materials provided with the distribution. 00014 * The names of its contributors may not be used to endorse or promote 00015 products derived from this software without specific prior written 00016 permission. 00017 00018 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 00019 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED 00020 TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A 00021 PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER 00022 OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 00023 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 00024 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 00025 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 00026 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 00027 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 00028 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00029 */ 00030 00031 #include "keypoint_interpolation.h" 00032 00033 namespace asrl { 00034 00035 00036 #define MID_IDX 1 00037 __global__ void fh_interp_extremum(float * d_hessian, Keypoint * d_features, int4 * d_maxmin, unsigned int * d_feature_counter, unsigned int * d_max_min_counter) 00038 { 00039 if(blockIdx.x >= *d_max_min_counter) 00040 return; 00041 00042 int octave = d_maxmin[blockIdx.x].w; 00043 00044 unsigned int i = ASRL_SURF_MAX_FEATURES; 00045 int hidx_x = d_maxmin[blockIdx.x].x-1 + threadIdx.x; 00046 int hidx_y = d_maxmin[blockIdx.x].y-1 + threadIdx.y; 00047 int hidx_z = d_maxmin[blockIdx.x].z-1 + threadIdx.z; 00048 00049 int hidx_lin = hidx_x + d_hessian_stride[0] * hidx_y + d_hessian_stride[0] * d_octave_params[octave].y_size * hidx_z; 00050 00051 __shared__ float fh_vals[3][3][3]; 00052 __shared__ Keypoint p; 00053 00054 fh_vals[threadIdx.z][threadIdx.y][threadIdx.x] = d_hessian[hidx_lin]; 00055 __syncthreads(); 00056 00057 if(threadIdx.x == 0 && threadIdx.y == 0 && threadIdx.z == 0){ 00058 __shared__ float H[3][3]; 00059 00060 //dxx 00061 H[0][0] = fh_vals[MID_IDX ][MID_IDX + 1][MID_IDX ] 00062 - 2.0f*fh_vals[MID_IDX ][MID_IDX ][MID_IDX ] 00063 + fh_vals[MID_IDX ][MID_IDX - 1][MID_IDX ]; 00064 00065 //dyy 00066 H[1][1] = fh_vals[MID_IDX ][MID_IDX ][MID_IDX + 1] 00067 - 2.0f*fh_vals[MID_IDX ][MID_IDX ][MID_IDX ] 00068 + fh_vals[MID_IDX ][MID_IDX ][MID_IDX - 1]; 00069 //dss 00070 H[2][2] = fh_vals[MID_IDX + 1][MID_IDX ][MID_IDX ] 00071 - 2.0f*fh_vals[MID_IDX ][MID_IDX ][MID_IDX ] 00072 + fh_vals[MID_IDX - 1][MID_IDX ][MID_IDX ]; 00073 00074 //dxy 00075 H[0][1]= 0.25f*(fh_vals[MID_IDX ][MID_IDX + 1][MID_IDX + 1] - 00076 fh_vals[MID_IDX ][MID_IDX - 1][MID_IDX + 1] - 00077 fh_vals[MID_IDX ][MID_IDX + 1][MID_IDX - 1] + 00078 fh_vals[MID_IDX ][MID_IDX - 1][MID_IDX - 1]); 00079 //dxs 00080 H[0][2]= 0.25f*(fh_vals[MID_IDX + 1][MID_IDX + 1][MID_IDX ] - 00081 fh_vals[MID_IDX + 1][MID_IDX - 1][MID_IDX ] - 00082 fh_vals[MID_IDX - 1][MID_IDX + 1][MID_IDX ] + 00083 fh_vals[MID_IDX - 1][MID_IDX - 1][MID_IDX ]); 00084 00085 //dys 00086 H[1][2]= 0.25f*(fh_vals[MID_IDX + 1][MID_IDX ][MID_IDX + 1] - 00087 fh_vals[MID_IDX + 1][MID_IDX ][MID_IDX - 1] - 00088 fh_vals[MID_IDX - 1][MID_IDX ][MID_IDX + 1] + 00089 fh_vals[MID_IDX - 1][MID_IDX ][MID_IDX - 1]); 00090 00091 //dyx = dxy 00092 H[1][0] = H[0][1]; 00093 00094 //dsx = dxs 00095 H[2][0] = H[0][2]; 00096 00097 //dsy = dys 00098 H[2][1] = H[1][2]; 00099 00100 00101 __shared__ float dD[3]; 00102 00103 //dx 00104 dD[0] = 0.5f*(fh_vals[MID_IDX ][MID_IDX + 1][MID_IDX ] - 00105 fh_vals[MID_IDX ][MID_IDX - 1][MID_IDX ]); 00106 //dy 00107 dD[1] = 0.5f*(fh_vals[MID_IDX ][MID_IDX ][MID_IDX + 1] - 00108 fh_vals[MID_IDX ][MID_IDX ][MID_IDX - 1]); 00109 //ds 00110 dD[2] = 0.5f*(fh_vals[MID_IDX + 1][MID_IDX ][MID_IDX ] - 00111 fh_vals[MID_IDX - 1][MID_IDX ][MID_IDX ]); 00112 00113 __shared__ float invdet; 00114 invdet = 1.f / 00115 ( 00116 H[0][0]*H[1][1]*H[2][2] 00117 + H[0][1]*H[1][2]*H[2][0] 00118 + H[0][2]*H[1][0]*H[2][1] 00119 - H[0][0]*H[1][2]*H[2][1] 00120 - H[0][1]*H[1][0]*H[2][2] 00121 - H[0][2]*H[1][1]*H[2][0] 00122 ); 00123 00124 00125 // // 1-based entries of a 3x3 inverse 00126 // /* [ |a22 a23| |a12 a13| |a12 a13|] */ 00127 // /* [ |a32 a33| -|a32 a33| |a22 a23|] */ 00128 // /* [ ] */ 00129 // /* [ |a21 a23| |a11 a13| |a11 a13|] */ 00130 // /* A^(-1) = [-|a31 a33| |a31 a33| -|a21 a23|] / d */ 00131 // /* [ ] */ 00132 // /* [ |a21 a22| |a11 a12| |a11 a12|] */ 00133 // /* [ |a31 a32| -|a31 a32| |a21 a22|] */ 00134 00135 __shared__ float Hinv[3][3]; 00136 Hinv[0][0] = invdet*(H[1][1]*H[2][2]-H[1][2]*H[2][1]); 00137 Hinv[0][1] = -invdet*(H[0][1]*H[2][2]-H[0][2]*H[2][1]); 00138 Hinv[0][2] = invdet*(H[0][1]*H[1][2]-H[0][2]*H[1][1]); 00139 00140 Hinv[1][0] = -invdet*(H[1][0]*H[2][2]-H[1][2]*H[2][0]); 00141 Hinv[1][1] = invdet*(H[0][0]*H[2][2]-H[0][2]*H[2][0]); 00142 Hinv[1][2] = -invdet*(H[0][0]*H[1][2]-H[0][2]*H[1][0]); 00143 00144 Hinv[2][0] = invdet*(H[1][0]*H[2][1]-H[1][1]*H[2][0]); 00145 Hinv[2][1] = -invdet*(H[0][0]*H[2][1]-H[0][1]*H[2][0]); 00146 Hinv[2][2] = invdet*(H[0][0]*H[1][1]-H[0][1]*H[1][0]); 00147 00148 __shared__ float x[3]; 00149 00150 x[0] = -(Hinv[0][0]*(dD[0]) + Hinv[0][1]*(dD[1]) + Hinv[0][2]*(dD[2])); 00151 x[1] = -(Hinv[1][0]*(dD[0]) + Hinv[1][1]*(dD[1]) + Hinv[1][2]*(dD[2])); 00152 x[2] = -(Hinv[2][0]*(dD[0]) + Hinv[2][1]*(dD[1]) + Hinv[2][2]*(dD[2])); 00153 00154 if(fabs(x[0])< 1.f && fabs(x[1])< 1.f && fabs(x[2])< 1.f) { // if the step is within the interpolation region, perform it 00155 00156 // Get a new feature index. 00157 i = atomicInc(d_feature_counter,(unsigned int) -1); 00158 00159 if(i < ASRL_SURF_MAX_FEATURES) { 00160 00161 //d_features[MAX_N_FEATURES*SF_X + i] = ((float)d_maxmin[blockIdx.x].x+x[1]) * (float)d_octave_params[octave].step + d_octave_params[octave].border; 00162 p.x = ((float)d_maxmin[blockIdx.x].x+x[1]) * (float)d_octave_params[octave].step + d_octave_params[octave].border; 00163 //d_features[MAX_N_FEATURES*SF_Y + i] = ((float)d_maxmin[blockIdx.x].y+x[0]) * (float)d_octave_params[octave].step + d_octave_params[octave].border; 00164 p.y = ((float)d_maxmin[blockIdx.x].y+x[0]) * (float)d_octave_params[octave].step + d_octave_params[octave].border; 00165 00166 if(x[2] > 0) 00167 { 00168 float a = d_hessian_scale[ASRL_SURF_MAX_INTERVALS*octave + d_maxmin[blockIdx.x].z]; 00169 float b = d_hessian_scale[ASRL_SURF_MAX_INTERVALS*octave + d_maxmin[blockIdx.x].z + 1]; 00170 //d_features[MAX_N_FEATURES*SF_SCALE+i] = (1.f - x[2]) * a + x[2] * b; 00171 p.size = (1.f - x[2]) * a + x[2] * b; 00172 } 00173 else 00174 { 00175 float a = d_hessian_scale[ASRL_SURF_MAX_INTERVALS*octave + d_maxmin[blockIdx.x].z]; 00176 float b = d_hessian_scale[ASRL_SURF_MAX_INTERVALS*octave + d_maxmin[blockIdx.x].z - 1]; 00177 //d_features[MAX_N_FEATURES*SF_SCALE+i] = (1.f + x[2]) * a - x[2] * b; 00178 p.size = (1.f + x[2]) * a - x[2] * b; 00179 } 00180 00181 //((char*)&d_features[MAX_N_FEATURES*SF_OCTAVE])[i] = octave; 00182 p.octave = octave; 00183 00184 //d_features[MAX_N_FEATURES*SF_STRENGTH + i] = fh_vals[MID_IDX][MID_IDX][MID_IDX]; 00185 p.response = fh_vals[MID_IDX][MID_IDX][MID_IDX]; 00186 00187 // Should we split up this transfer over many threads? 00188 d_features[i] = p; 00189 } else { 00190 *d_feature_counter = ASRL_SURF_MAX_FEATURES; 00191 } // If there is still room in memory for the feature 00192 } // If the subpixel interpolation worked 00193 } // If this is thread 0. 00194 00195 // __syncthreads(); 00196 // // Transfer the feature back. 00197 // if(i < MAX_N_FEATURES) 00198 // { 00199 // int tid = threadIdx.x + __mul24(threadIdx.y,3) + __mul24(threadIdx.z,9); 00200 // if(tid < 6) 00201 // { 00202 // ((float*)&d_features[i])[tid] = ((float*)&p)[tid]; 00203 // } 00204 // } 00205 } 00206 00207 00208 void run_fh_interp_extremum(float * d_hessian, Keypoint * d_features, int4 * d_maxmin, unsigned int * d_feature_counter, unsigned int * d_max_min_counter) 00209 { 00210 dim3 threads; 00211 threads.x = 3; 00212 threads.y = 3; 00213 threads.z = 3; 00214 00215 dim3 grid; 00216 grid.x = ASRL_SURF_MAX_CANDIDATES; 00217 grid.y = 1; 00218 grid.z = 1; 00219 00220 fh_interp_extremum<<< grid, threads >>>(d_hessian, 00221 d_features, 00222 d_maxmin, 00223 d_feature_counter, 00224 d_max_min_counter); 00225 } 00226 00227 } // namespace asrl
Generated on Fri Apr 30 20:06:20 2010 for gpusurf by
1.6.2