Random walk within ROI of an image

Question

I'm working on a diffusion MRI project and have a c++ problem which I need help with.

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.

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)

Code of currently used method:

blitz::Array<Catmull,1> FIBERs;         // Fiber with control nodes (knots) which I want to move
NIFTI<INT16>* niiWM = new NIFTI<INT16>; // .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 >= 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<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->img)(Vx,Vy,Vz) == 0){              // Penalize when outside the WM
    Pe += 1;}
     }
     return Pe;     
}

I hope that this was enough information. Thanks for your time.

Answer 1

What type of cost function do you have (if I am reading from your code it is just the summation of pixels which do not fall outside the white matter?)? If it is smooth and continuous, you could try gradient descent to find the minimum (since yours looks discrete, though, is there a way for you to approximate gradient descent?). This should function more efficiently than a random walk.

If maintaining randomness is important, then the best way to do that without suffering a lot of re-draws is to re-normalize the range of distribution you can select from. If I understand correctly, you have the white matter mask already, so one method that should allow you to do this would be:

1.) Assign a probability to all white matter pixels based on a Gaussian centered at your current location

2.) Re-normalize those probabilities so they sum to 1

3.) Pull a random number from a uniform distribution and find which pixel it corresponds to

4.) Move your control point there

So, for example, if you have four possible white matter pixels you could move to with probability 0.1, 0.17, 0.2, and 0.23, at respective locations (x1,y1), (x2,y2), (x3, y3), (x4, y4) with all other probable pixels (based on your Gaussian random step) being non-white matter pixels that will force you to re-draw, then you would re-scale your probabilities by dividing by their sum (0.7 in this case). Then pull a random number from a uniform distribution (let's say you pulled 0.3). This would correspond to the second point (x2, y2) based on the re-normalized probabilities, and you would move your control point there.

As you can see, this has a fair number of computational steps to perform, but will ensure you only pull one random target to update; any improvement in performance will therefore depend on the relative cost of these extra computations versus how often you are actually re-sampling. I also wrote this up pretty fast, so you might be able to come up with a more efficient way to limit the scope of your random selection to only the good pixels.