Humboldt-Universität zu Berlin - Faculty of Mathematics and Natural Sciences - Strukturforschung / Elektronenmikroskopie

SelectAndAlign2Images.txt

// This script finds a sub-region of one image within another image:

Image img1, img2, chi2Map, imgSmall, imgBig, imgShifted
number width,height,top=0,left=0,bottom,right,w1,h1,w2,h2,top2,bottom2,left2,right2
ROI subRegion

///////////////////////////////////////////////////////////
// This function finds the position of a small image within
// a bigger one:
///////////////////////////////////////////////////////////
// The images are passed by reference (images always are):
image findSmallInBig(image &big,image &small) {

number Nx,Ny,Nx2,Ny2,Nx3,Ny3;
complexImage chi2Map,timg1, mask, subMask;
Image realChi2Map

GetSize(big,Nx,Ny);
GetSize(small,Nx2,Ny2);

// allocate memory for the temporary image:
timg1 := ComplexImage("reciprocal space Image", 8, Nx, Ny)

// create a mask that has the same size as the small image:
mask = ExprSize(Nx,Ny,0)
mask[0,0,Ny2,Nx2] = 1


// compute the FFT of mask:
T_fft_C2C(mask,mask)
// compute the FFT of big^2 and save it in timg1:
timg1 = big*big;
T_fft_C2C(timg1,timg1)

chi2Map = conjugate(mask)*timg1

mask[] = 0;
// assign a variable name to the sub-region of the mask that we want to copy
// the small image to:
mask[0,0,Ny2,Nx2] = small


T_fft_C2C(mask,mask);
timg1 = big
T_fft_C2C(timg1,timg1);
chi2Map = chi2Map-2*conjugate(mask)*timg1;
T_ifft_C2C(chi2Map,chi2Map);

// shift the center of the chi2-map to the center of the image:
chi2Map := T_shiftImageCenterComplex(chi2Map)
// scale and convert the chi2Map to a real image:
// realChi2Map = real(chi2Map)/(Nx*Ny) + sum(small*small)
realChi2Map = real(chi2Map)
// showimage(realChi2Map)

// free the allocated space again:
deleteImage(chi2Map)
deleteImage(timg1)
deleteImage(mask)

return realChi2Map;
}
///////////////////////////////////////////////////////////


// Ask the user, which 2 images he want to have correlated:
GetTwoImagesWithPrompt("select 2 images:","Image correlation", img1,img2)

// Find out, how big these images are:
img1.GetSize(w1,h1)
img2.GetSize(w2,h2)


if (w1*h1 != w2*h2) {
// if one image is smaller than the other one, we will try to find the
// smaller one within the bigger image
if (w1*h1 < w2*h2) {
imgSmall:=img1
imgBig := img2
imgShifted := img1
}
else {
imgSmall:=img2
imgBig := img1
imgShifted := img2
}
// Just in case the images have totally different shapes:
if (w1 > w2) imgSmall := imgSmall[0,0,h1,w2];
if (h1 > h2) imgSmall := imgSmall[0,0,h2,w1];
}
else {
// If both images have the same size, we will try to find a selection
// in one of the images
ImageDisplay imgDisp = img1.ImageGetImageDisplay(0)
number ROIcount = ImageDisplayCountROIs(imgDisp);

if (ROIcount == 0) {
imgDisp = img2.ImageGetImageDisplay(0)
ROIcount = ImageDisplayCountROIs(imgDisp);
if (ROIcount == 0) {
GetSize(img2,width,height)
subRegion = CreateROI();
Result("No region of interest has been selected - will use central region of image 2\n")
ROISetRectangle(subRegion,floor(0.25*height),floor(0.25*width),ceil(0.75*height),ceil(0.75*width));
ImageDisplayAddROI(imgDisp,subRegion);
}
subRegion = imgDisp.ImageDisplayGetROI(0)
subRegion.ROIGetRectangle(top, left, bottom, right)
imgSmall := img2[top,left,bottom,right]
imgBig := img1
imgShifted := img2
}
else {
subRegion = imgDisp.ImageDisplayGetROI(0)
subRegion.ROIGetRectangle(top, left, bottom, right)
imgSmall = img1[top,left,bottom,right]
imgBig := img2
imgShifted := img1
}
}
chi2Map=findSmallInBig(imgBig,imgSmall)

// find the relative image shift:
number mx,my
number mval = min(chi2Map,mx,my)
chi2Map.getSize(width,height)

mx = mx-floor(width/2);
my = my-floor(height/2);

result("Position of minimum in chi2-map: "+mx+", "+my+" (val="+mval+")\n");

///////////////////////////////////////////////////////////////////////
// We must now add the 2 images according to the relative image shift:

imgShifted.getSize(w1,h1);
imgBig.getSize(w2,h2);

// In order to add the 2 images, we must make sure that we keep track of where
// the edges are
top2 = min(-my-top,0)
left2 = min(-mx-left,0)
bottom2 = max(h1+abs(my+top),h2)
right2 = max(w1+abs(mx+left),w2)
// The sum image must have the maximum size:
Image imgSum := RealImage("Sum image",8,right2,bottom2)

// add both image to the sum image
imgSum[-top2,-left2,h2-top2,w2-left2] += imgBig
imgSum[-top2-my-top,-left2-mx-left,-top2-my+h1-top,-left2-mx+w1-left] += imgShifted
showimage(imgSum)