title: Improving C++ algorithm for finding all points within a sphere of radius r
body: <p>Language/Compiler: C++ (Visual Studio 2013)
  Experience: ~2 months</p>
  
  <p>I am working in a rectangular grid in 3D-space (size: xdim by ydim by zdim) where , "xgrid, ygrid, and zgrid" are 3D arrays of the x,y, and z-coordinates, respectively. Now, I am interested in finding all points that lie within a sphere of radius "r" centered about the point "(vi,vj,vk)". I want to store the index locations of these points in the vectors "xidx,yidx,zidx". For a single point this algorithm works and is fast enough but when I wish to iterate over many points within the 3D-space I run into very long run times. </p>
  
  <p>Does anyone have any suggestions on how I can improve the implementation of this algorithm in C++? After running some profiling software I found online (very sleepy, Luke stackwalker) it seems that the "std::vector::size" and "std::vector::operator[]" member functions are bogging down my code. Any help is greatly appreciated.</p>
  
  <p>Note: Since I do not know a priori how many voxels are within the sphere, I set the length of vectors xidx,yidx,zidx to be larger than necessary and then erase all the excess elements at the end of the function.</p>
  
  <pre><code>void find_nv(int vi, int vj, int  vk, vector&lt;double&gt; &amp;xidx, vector&lt;double&gt; &amp;yidx, vector&lt;double&gt; &amp;zidx, double*** &amp;xgrid, double*** &amp;ygrid, double*** &amp;zgrid, int r, double xdim,double ydim,double zdim, double pdim)
  
  {
  double xcor, ycor, zcor,xval,yval,zval;
  vector&lt;double&gt;xyz(3);
  xyz[0] = xgrid[vi][vj][vk];
  xyz[1] = ygrid[vi][vj][vk];
  xyz[2] = zgrid[vi][vj][vk];
  int counter = 0;
  
  // Confine loop to be within boundaries of sphere
  int istart = vi - r;
  int iend = vi + r;
  int jstart = vj - r;
  int jend = vj + r;
  int kstart = vk - r;
  int kend = vk + r;
  
  if (istart &lt; 0) {
      istart = 0;
  }
  if (iend &gt; xdim-1) {
      iend = xdim-1;
  }
  if (jstart &lt; 0) {
      jstart = 0;
  }
  if (jend &gt; ydim - 1) {
      jend = ydim-1;
  }
  if (kstart &lt; 0) {
      kstart = 0;
  }
  if (kend &gt; zdim - 1)
      kend = zdim - 1;
  
  //-----------------------------------------------------------
   // Begin iterating through all points
  //-----------------------------------------------------------
  for (int k = 0; k &lt; kend+1; ++k)
  {
      for (int j = 0; j &lt; jend+1; ++j)
      {
          for (int i = 0; i &lt; iend+1; ++i)
          {
              if (i == vi &amp;&amp; j == vj &amp;&amp; k == vk)
                  continue;
              else
              {
              xcor = pow((xgrid[i][j][k] - xyz[0]), 2);
              ycor = pow((ygrid[i][j][k] - xyz[1]), 2);
              zcor = pow((zgrid[i][j][k] - xyz[2]), 2);
              double rsqr = pow(r, 2);
              double sphere = xcor + ycor + zcor;
              if (sphere &lt;= rsqr)
              {
                  xidx[counter]=i;
                  yidx[counter]=j;
                  zidx[counter] = k;
                  counter = counter + 1;
              }
  
              else
              {
              }
              //cout &lt;&lt; "counter = " &lt;&lt; counter - 1;
              }
          }
      }
  }
  
  // erase all appending zeros that are not voxels within sphere
  xidx.erase(xidx.begin() + (counter), xidx.end()); 
  yidx.erase(yidx.begin() + (counter), yidx.end()); 
  zidx.erase(zidx.begin() + (counter), zidx.end()); 
  
  return 0;
  </code></pre>