-1

I have two 3D images, i need to register these two images using "lsqcurvefit". I know that I can use "imregister" but I want to use my own registration using "lsqcurvefit" in Matlab. My images are are following Gaussian distribution. it is not documented well that how should I provide it, anyone can help me in detail?

image registration is a repeated process of maping source image to target image using i.e affine. i want to use intensity base registration, and i use all voxels of my image. therefore, i need to fit these two images as much as possible.

Thanks

Ehsan
  • 517
  • 1
  • 7
  • 32
  • As in your image is a single Gaussian peak? Or multiple Gaussian peaks? – Staus Feb 24 '15 at 05:08
  • @Staus multiple peaks – Ehsan Feb 24 '15 at 05:11
  • 1
    Quick version - use something like `imregionalmax' to find local maxima, take a sub-image around each of a defined size (like expected sigma*3 or so), then fit each to a 2d Gaussian using lsqcurvefit. Given the resulting list of localizations in each image you can then use lsqcurvefit to fit the transform matrix ([cos theta, -sin theta, translate x; sin theta, cos theta, translate y; 0 0 1]) in the equation LocalizationsB = TransformMatrix * Localizations A. – Staus Feb 24 '15 at 05:32
  • @Staus thank for your reply, what if I have Single gaussian peak? i found out that some of my images are single gaussian peak , thanks in advance – Ehsan Feb 24 '15 at 08:34
  • 1
    You need at least two points to define a transform with rotation and translation. Given one you can do translation only, so change the matrix above to [1 0 TransX; 0 1 TransY; 0 0 1] in the cases where you have one point, or go for a more image-based registration. – Staus Feb 24 '15 at 19:29
  • @Staus I could fit the result of local maxima with 2d Gaussian using lsqcurvefit. but i dont know hot to fit transform matrix. can you help me in this part? Thanks – Ehsan Feb 25 '15 at 01:26

1 Answers1

1

Here's an example of how to do point-wise image registration using lsqcurvefit. Basically you make a function that takes a set of points and an Affine matrix (we're just going to use the translate and rotate parts but you can use skew and magnify if desired) and returns a new set of points. There's probably a built-in function for this already but it's only two lines so it's easy to write. That function is:

function TformPts = TransformPoints(StartCoordinates, TransformMatrix)

TformPts = StartCoordinates*TransformMatrix;

Here's a script that generates some points, rotates and translates them by a random angle and vector, then uses the TransformPoints function as the input for lsqcurvefit to fit the needed transformation matrix for the registration. Then it's just a matrix multiplication to generate the registered set of points. If we did this all right the red circles (original data) will line up with the black stars (shifted then registered points) very well when the code below is run.

% 20 random points in x and y between 0 and 100
% row of ones pads out third dimension
pointsList = [100*rand(2, 20); ones(1, 20)];

rotateTheta = pi*rand(1); % add rotation, in radians
translateVector = 10*rand(1,2); % add translation, up to 10 units here

% 2D transformation matrix
% last row pads out third dimension
inputTransMatrix = [cos(rotateTheta), -sin(rotateTheta), translateVector(1);
                    sin(rotateTheta), cos(rotateTheta), translateVector(2);
                    0 0 1];


% Transform starting points by this matrix to make an array of shifted
% points.  
% For point-wise registration, pointsList represents points from one image,
% shiftedPoints points from the other image
shiftedPoints = inputTransMatrix*pointsList;
% Add some random noise
% Remove this line if you want the registration to be exact
shiftedPoints = shiftedPoints + rand(size(shiftedPoints, 1), size(shiftedPoints, 2));

% Plot starting sets of points
figure(1)
plot(pointsList(1,:), pointsList(2,:), 'ro');
hold on
plot(shiftedPoints(1,:), shiftedPoints(2,:), 'bx');
hold off

% Fitting routine
% Make some initial, random guesses
initialFitTheta = pi*rand(1);
initialFitTranslate = [2, 2];

guessTransMatrix = [cos(initialFitTheta), -sin(initialFitTheta), initialFitTranslate(1);
                    sin(initialFitTheta), cos(initialFitTheta), initialFitTranslate(2);
                    0 0 1];

% fit = lsqcurvefit(@fcn, initialGuess, shiftedPoints, referencePoints)             
fitTransMatrix = lsqcurvefit(@TransformPoints, guessTransMatrix, pointsList, shiftedPoints);

% Un-shift second set of points by fit values
fitShiftPoints = fitTransMatrix\shiftedPoints;

% Plot it up
figure(1)
hold on
plot(fitShiftPoints(1,:), fitShiftPoints(2,:), 'k*');
hold off

% Display start transformation and result fit
disp(inputTransMatrix)
disp(fitTransMatrix)
Staus
  • 651
  • 3
  • 7