title: Random walk within ROI of an image
body: <p>I'm working on a diffusion MRI project and have a c++ problem which I need help with. </p>
  
  <p>Background:
  I have a .nii image as a reference in which I would like a control node (which is part of a spline) to perform a random walk within a constrained ROI of the image. To clarify: the .nii image is a mask showing the brains white matter (WM), I want a spline with 4 control nodes to randomly move within this mask until they reach a state which minimizes a cost function. Currently what I'm doing is just updating a random control node to move in a random direction according to a Gaussian distribution and then checking if it is inside or outside the WM and redrawing if is outside the WM. But this can lead to many unnecessary draws which I would like to avoid.</p>
  
  <p>To my question: 
  Is there any better way to do this? I have an idea about collecting all the pixel coordinates of my ROI into a separate array or matrix which I can draw from the indexes of that matrix instead? This method would also be helpful if I want to reinitialize my splines with four completely new control points. (Which I also need to implement in the future)</p>
  
  <p>Code of currently used method:</p>
  
  <pre><code>blitz::Array&lt;Catmull,1&gt; FIBERs;         // Fiber with control nodes (knots) which I want to move
  NIFTI&lt;INT16&gt;* niiWM = new NIFTI&lt;INT16&gt;; // .nii Image which contains WM mask
  void SimulatedAnnealing_OneStep( void )
  {
      iF = floor( FIBERs.extent(0)*uniformGen.random());  //  select randomly a FIBER 
      NewProposal(iF);
      Pen = CheckWM(iF);
      while(Pen &gt;= 1){
          NewProposal(iF);
      Pen = CheckWM(iF);
      } 
  }
  void NewProposal( int ifff)
  {
      static int              iK;
      static POINT            delta;
  
      delta.Set( MOVE_sigma*normalGen.random(), MOVE_sigma*normalGen.random(), 0 );       
      iK = floor(uniformGen.random()*6); 
  
      if( iK==2 || iK==3 ) // Moving the middle control points within the WM 
      {   // FIBERs(ifff) is a randomly selected fiber (outside this method) which we are moving
          FIBERs(ifff).KNOTs[iK].x += delta.x;
          FIBERs(ifff).KNOTs[iK].y += delta.y;
          FIBERs(ifff).KNOTs[iK].z += delta.z;
          FIBERs(ifff).KNOTs[iK].x =  fmin( fmax( FIBERs(ifff).KNOTs[iK].x,0), Nx );
          FIBERs(ifff).KNOTs[iK].y =  fmin( fmax( FIBERs(ifff).KNOTs[iK].y,  0), Ny );
          FIBERs(ifff).KNOTs[iK].z =  fmin( fmax( FIBERs(ifff).KNOTs[iK].z, 0), Nz ); 
      }   
  }
  
  int CheckWM( int iff )
  {
  int j,Pe;
  int Vx, Vy, Vz;
  
  Pe = 0;
  
  for(j=0; j&lt;FIBERs(iff).P.extent(0) ;j++){           // Loop through the segments on the curve
      Vx = floor( FIBERs(iff).P(j).x );               // midpoint of the segment
      Vy = floor( FIBERs(iff).P(j).y );
      Vz = floor( FIBERs(iff).P(j).z );
  
      if ((*niiWM-&gt;img)(Vx,Vy,Vz) == 0){              // Penalize when outside the WM
      Pe += 1;}
       }
       return Pe;     
  }
  </code></pre>
  
  <p>I hope that this was enough information. Thanks for your time.</p>