Héctor here, your AFM expert at Nanosurf calling out for people to share their Friday afternoon experiments. Today I will show you how to combine multiple scans to obtain an image larger than the one accesible by the maximum scanner range. (Mini tutorial of MountainsSPIP and Python).
On many fridayAFM posts I showed you images that are a combination of multiple individual scans. Examples of that are for instance the white rose petal:
the pine seed:
the 1982 floppy disc:
and of course, the fridayAFM - 1984 EPROM:
All these have been made with MountainsSPIP, using the Patch tool and manual alignment. The process is relatively simple, and I'll show you how to do it in a moment, however, there are several cases where it is almost impossible to achieve it (this is why combining it with Python, and the reason of this fridayAFM).
In order to continue my study of magnetic media bit size vs production year (see more in this Nanosurf Insights post), I obtained a 1973 8 inch floppy disc. The problem here is that the track size is more than 100 μm, and there are no visible features to align the individual scans (at least I am not able to distinguish them).
So, lets see how to stitch images using MountainsSPIP (I will use the rose petal as example because it is easier to align), then we will see how to help us with Python, and how by combining both we can align the floppy disc images and obtain the bit size.
1. - Load the data onto MountainsSPIP. Simply by opening the software and dragging the file onto the working area.
This will automatically generate a visualization of the data, and it works for most of the AFM data formats outthere.
You can change the size of this visualization and select other channels a part from the topography. It is also possible to change visual elements like displaying the xis or not, the background... or the colormap (always go for Gold if possible, people are more used to understand features in that colormap).
In my case, I don't like this visualization a lot and I prefer having the different channels separated, so I right click on the workflow element corresponding to the visualization, and delete it (don't worry, this doesn't delete the data).
Next, with the data selected in the workflow, I use the Extract channels operator.
In this case only the topography is interesting, but in the case of the floppy disc I will extract topography and phase signal during the second pass.
The result looks something like this:
2. - Data correction. Once the data has been loaded onto the software, it is time to process it. Using the workflow, select the step in the flow where do you want to insert an operator. In this case, select the pseudo-color view and insert a level line by line (one of the most basic corrections).
This will correct offset shifts between scan lines or the tilt of the sample, however, in this case, the data seems good enough as it is and the process varely produces any change in the visualization.
Note that you can use the side menus to change the color scale and improve the visualization of defects. Once this step is completed, it will appear in the workflow. There you can choose to add more steps before or after, or to create more visualizations from the same data.
It is important to notice that although you can change the color scale directly on the lateral bar of the pseudo-color plot:
However, the best approach to keep the scale linear and to control the limits is to go to the axis menu.
I usually go for simple numbers without many decimals.
3. - Exporting. Once all the data manipulation is completed (e.g. removing the axes and the color scale if one wants to further manipulate the data in Python), there are several options to export the data. The most basic ones are to save the document as a pdf.
or to export the pseudo-color view alone as an image.
4. - Data analysis. This is not strictly necessary for image stitching, but could be useful as part of this introduction to MountainsSPIP. One of the most basic data analysis process is extracting cross section. This can be done by selecting the pseudo-color view and choosing extragin a profile from the Operators menu.
Although the profile extraction menu is quite complete, my recommendation, to ease interpretation and analysis, is to restrict ourselves to straight lines. It is also sometimes recommended to take several profiles and average them to reduce the noise (that is done by changing the half width parameter).
One nice feature is that these profiles are dynamic, meaning that you can move the line around and see how the profile changes
The next type of analysis is roughness measurements. To do that, choose the leveled data and from the right click menu, choose parameters table.
This will open a menu that will let you choose from a large list of roughness indicators.
Last two tricks are that you can drag and drop another file instead of the original one and all the operation will be automatically applied to it (another way of achieving this will be selecting all the steps and in the minidocs tab click create a minidoc, this will create a custom analisis which consists in all the steps we just performed and that is ready to be applied to any new data you load onto the software).
The last trick is that in order to add text you just double click on an empty space and start typing (obvious, but sometimes hard to realize how it works).
By the way, notice (above) how this new data file needs some fine-tuning of the leveling, you can do that by right clicking on the operator in the workflow and clicking on Recall the operator.
5. - Image stitching. Once you have uploaded all the data and is nicely processed (you need to judge during the stitching if the data has to be leveled or not), the next step is to apply the Patch operator to stitch the images together. Why not the stitching operator? Stitching means no overlap, while patch means some overlap. In our case we will use the overlap to properly align the images, and in general, it is very likely (and recommended), that your images will overlap a bit, hence why the patch and not the stitch operator. By the way, before entering the stitching menu, it is important to write down the numbers that appear in the workflow next to the last step performed on each file. These numbers will be used to identify the images that will be combined together.
The idea behind the stitching menu is that you add images from the workflow by copying them, and the you can move each image in x,y, and z and you look at the result on the right visualization window.
The nice thing is that it can compute the z offset from the overlapping regions, making the z alignment easier, the bad thing is that the visualization window is very small and it becomes tedious to navigate around in it.
Usually the process I follow is this, I left click on the new image in the result window and move it around by holding the left mouse button. This process allows me to do a rough alignment, and then, using the arrows on the x and y offset I do a fine-tune while I zoomed in a small region on the result window. At this point is when you will notice if the leveling was good or if it should be removed completely. By iterating this process, I can align any number of images, but have in mind that the process is never going to be perfect, so there is a trade of between what you see and how many imperfections the resulting image has.
Note that it is possible to use the topography channel to align the data and the copy the x and y coordinates and align the phase signal (this is how I did the 1982 floppy disc), however, in the case of the 8-inch floppy disc, the topography doesn't have many recognizable features and the alignment process is very difficult.
So... this is where I thought, what if I use Python to obtain the x and y coordinates and then do the stitching in Mountains using that information? In principle I could also do the stitching in Python... but handling the 3D data alongside the alignment is almost like building another MountainsSPIP in Python (which I'm not interested in doing). So lets restrict ourselves to obtain the x and y coordinates in Python and leave the heavy lifting to MountainsSPIP.
For this next section we will try OpenCV, which is a Python library dedicated to computer vision.
I will start from the example shown in this tutorial (or this one) and modify it slightly for our needs.
Basically, the tutorial looks for key features in two images and then calculates the transformation needed on one of the images so the matched features overlap. This means, rotation, scaling, warping... Unfortunately, this is not suitable for our data, as one important thing is that scales and ratios are conserved. For this reason, my main modification to the script consists in replacing "find homography", which calculates the transformation needed, for "estimateAffine2D", which also calculates the transformation but only applying translations and rotations. This also means that we need to use "warpAffine" instead of "warpPerspective". Just in case you want to play with general stitching, I left the commands in there and you just need to comment/uncomment the mentioned commands. (If you want to play even more, here is an example of how to use homography to replace a billboard within an image with a custom image).
I also added another variation. The original code was meant for panorama formation, assuming the final image is cropped to the smallest combined field of view... In our case we don't want that, so I added an extra space around one of the images before stitching the second one, and it is in the end when I crop down the excess.
In order to find the offset in meters between images, I assume both images have the same pixel size (they can have different number of pixels as long as their pixels have the same size), I ask the user for the lateral size of one of the images... and then use a box fitting tool to identify the location of each image. I'm sure there are more elegant methods, which are less prone to error... but the code to achieve this is almost the same I used for cropping, so it was convenient.
Here is my code:
#%% # Import the necessary packages import numpy as np import imutils import cv2 from pathlib import Path import tkinter as tk from tkinter import filedialog # I personally don't use classes that much, but I don't want to modify a lot the example, so most of what appears here is copied from: # https://medium.com/analytics-vidhya/panorama-formation-using-image-stitching-using-opencv-1068a0e8e47b class Stitcher: def __init__(self): # determine if we are using OpenCV v3.X self.isv3 = imutils.is_cv3(or_better=True) def stitch(self, images, ratio=0.85, reprojThresh=5.0, showMatches=False): # unpack the images, then detect keypoints and extract # local invariant descriptors from them (imageB, imageA) = images ####################New################### # The original script proceeds with imageA and imageB, but if the warping goes left or top, the result is cropped, that is why I create imagebb which is basically B # but with a large empty border around as wide as imageA imagebb=np.zeros((imageB.shape[0]+2*imageA.shape[0],imageB.shape[1]+2*imageA.shape[1],3),dtype=np.uint8) imagebb[imageA.shape[0]:imageA.shape[0]+imageB.shape[0], imageA.shape[1]:imageA.shape[1]+imageB.shape[1],0] = imageB[:,:,0] imagebb[imageA.shape[0]:imageA.shape[0]+imageB.shape[0], imageA.shape[1]:imageA.shape[1]+imageB.shape[1],1] = imageB[:,:,1] imagebb[imageA.shape[0]:imageA.shape[0]+imageB.shape[0], imageA.shape[1]:imageA.shape[1]+imageB.shape[1],2] = imageB[:,:,2] ####################end of New################### (kpsA, featuresA) = self.detectAndDescribe(imageA) (kpsB, featuresB) = self.detectAndDescribe(imagebb) # match features between the two images M = self.matchKeypoints(kpsA, kpsB, featuresA, featuresB, ratio, reprojThresh) # if the match is None, then there aren't enough matched # keypoints to create a panorama if M is None: return None # otherwise, apply a perspective warp to stitch the images # together (matches, H, status) = M ####################New################### # If you want a general transformation, uncomment warpPerspective and comment warpAffine. # There is a piece of code that also needs commenting later on. # Perspective scales and warps the images to match the matched points, Affine only does translations and rotations. #result = cv2.warpPerspective(imageA, H, # (2*imageA.shape[1] + imageB.shape[1], 2*imageA.shape[0]+imageB.shape[0])) #print(H) result = cv2.warpAffine(imageA, H, (2*imageA.shape[1] + imageB.shape[1], 2*imageA.shape[0]+imageB.shape[0])) ####################end of New################### ####################New################### # At this point, ImageA has been moved to the final position (on a matrix as large as bb), # So we use the threshold and bonding box trick to find its position (I'm sure somebody more inteligent than me can figure out from H... if you are, please share back how to.) grayscale = cv2.cvtColor(result, cv2.COLOR_BGR2GRAY) th,thresholded = cv2.threshold(grayscale, 0, 255, cv2.THRESH_OTSU) bbox = cv2.boundingRect(thresholded) x, y, w, h = bbox distx=-(result.shape[1]/2-(x+w/2)) disty=(result.shape[0]/2-(y+h/2)) #The next step is to copy imageB onto the result in the middel of the array. result[imageA.shape[0]:imageA.shape[0]+imageB.shape[0], imageA.shape[1]:imageA.shape[1]+imageB.shape[1]] = imageB ####################end of New################### # check to see if the keypoint matches should be visualized if showMatches: vis = self.drawMatches(imageA, imagebb, kpsA, kpsB, matches, status) # return a tuple of the stitched image and the # visualization return (result, vis,distx, disty) # return the stitched image return (result, distx, disty) def detectAndDescribe(self, image): # convert the image to grayscale gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY) # check to see if we are using OpenCV 3.X if self.isv3: # detect and extract features from the image descriptor = cv2.xfeatures2d.SIFT_create() (kps, features) = descriptor.detectAndCompute(image, None) # otherwise, we are using OpenCV 2.4.X else: # detect keypoints in the image detector = cv2.FeatureDetector_create("SIFT") kps = detector.detect(gray) # extract features from the image extractor = cv2.DescriptorExtractor_create("SIFT") (kps, features) = extractor.compute(gray, kps) # convert the keypoints from KeyPoint objects to NumPy # arrays kps = np.float32([kp.pt for kp in kps]) # return a tuple of keypoints and features return (kps, features) def matchKeypoints(self, kpsA, kpsB, featuresA, featuresB, ratio, reprojThresh): # compute the raw matches and initialize the list of actual # matches matcher = cv2.DescriptorMatcher_create("BruteForce") rawMatches = matcher.knnMatch(featuresA, featuresB, 2) matches = [] # loop over the raw matches for m in rawMatches: # ensure the distance is within a certain ratio of each # other (i.e. Lowe's ratio test) if len(m) == 2 and m[0].distance < m[1].distance * ratio: matches.append((m[0].trainIdx, m[0].queryIdx)) # computing a homography requires at least 8 matches if len(matches) > 8: # construct the two sets of points ptsA = np.float32([kpsA[i] for (_, i) in matches]) ptsB = np.float32([kpsB[i] for (i, _) in matches]) ####################New################### #Either computes the homography or the Affine translation, as mentioned above. Comment/uncomment depending which one you want to use # Compute the homography between the two sets of points #(H, status) = cv2.findHomography(ptsA, ptsB, cv2.RANSAC, # reprojThresh) # Instead calculates a rigid transformation H, status = cv2.estimateAffine2D(ptsA, ptsB) ####################end of New################### # return the matches along with the homograpy matrix # and status of each matched point return (matches, H, status) # otherwise, no homograpy could be computed return None def drawMatches(self, imageA, imageB, kpsA, kpsB, matches, status): # initialize the output visualization image (hA, wA) = imageA.shape[:2] (hB, wB) = imageB.shape[:2] vis = np.zeros((max(hA, hB), wA + wB, 3), dtype="uint8") vis[0:hA, 0:wA] = imageA vis[0:hB, wA:] = imageB # loop over the matches for ((trainIdx, queryIdx), s) in zip(matches, status): # only process the match if the keypoint was successfully # matched if s == 1: # draw the match ptA = (int(kpsA[queryIdx][0]), int(kpsA[queryIdx][1])) ptB = (int(kpsB[trainIdx][0]) + wA, int(kpsB[trainIdx][1])) cv2.line(vis, ptA, ptB, (0, 255, 0), 1) # return the visualization return vis #GUI to select the datafile root = tk.Tk() root.withdraw() #Assumes images have the same scale. print('Input Image A') filenameA = filedialog.askopenfilename() name=Path(filenameA).stem folder=directory = os.path.split(filenameA) imageA = cv2.imread(filenameA) print('Input Image B') filenameB = filedialog.askopenfilename() imageB = cv2.imread(filenameB) # If the wresolution is too high the computer can be very slow or fail completely, in my case, if the images have more than 6000 pixels, I scale them down (hA, wA) = imageA.shape[:2] if wA>6000: imageA = imutils.resize(imageA, width=int(wA/10)) (hB, wB) = imageB.shape[:2] imageB = imutils.resize(imageB, width=int(wB/10)) # This is to calculate the offset in um widthA=input ('Image A width in um') singlepixelwidth=float(widthA)/imageA.shape[1] stitcher = Stitcher() ####### This is the main process (result, vis,distx,disty) = stitcher.stitch([imageA, imageB], showMatches=True) ####### print('Relative x movement in um of image B', distx*singlepixelwidth) print('Relative y movement in um of image B', disty*singlepixelwidth) # Display the images # First the input images cv2.imshow("Image A", imageA) cv2.imshow("Image B", imageB) # Then the keypoints matched cv2.imshow("Keypoint Matches", vis) # Then we need to crop the result to eliminate the part of the image needed for the process but not wanted for the final result. # Convert to gray grayscale = cv2.cvtColor(result, cv2.COLOR_BGR2GRAY) cv2.imshow("Grayscale",grayscale) # Apply a threshold to identify the empty parts th,thresholded = cv2.threshold(grayscale, 0, 255, cv2.THRESH_OTSU) cv2.imshow("Thresholded",thresholded) # Create the smallest bounding box that contains all the data in the image bbox = cv2.boundingRect(thresholded) x, y, w, h = bbox # The final data, called foreground is the area defined by the bounding box. foreground = result[y:y+h, x:x+w] cv2.imwrite(str(folder)+'/save_image_name.png',foreground) cv2.imshow("Result", result) cv2.imshow("Result_Cropped", foreground) cv2.waitKey(0) #%%
The outputs of the code are the two input images. In this case two images of the rose petal.
The image showing the matches found between the images.
The result of stitching both images.
The grayscale and the threshold image (this is just to check if the cropping of the image is working properly).
Here the thresholded image:
and finally, the final result alongside with the offset from one image to another, which are -68 μm in x and 1.8 μm in y.
These are the numbers we have to input in MountaisSPIP, and if we do, we find that the match works.
You will have to write them down and add together successive displacements, when combining more than one image, but as far as I can tell the numbers are good enough to not require fine tuning by hand.
So... now, can we use this process with the topography images of the 1973 floppy disc and use the resulting numbers to stitch the phase images? Judge by yourselves.
This gives us a bit size of 22.5 by 289 μm2 . Which once added to the bit size over the years graph...
... opens new questions. Is it an outliner? Although the technology was introduced in 1973, did the bit size changed over the years? Or was the curve flat until in the 80s an accelerated growth happened? Maybe we will be able to answer these questions in future fridayAFMs.
Let's recap. We have seen how to stitch images using MountainsSPIP, alongside some basic operations like flattening, profile extraction or roughness measurements. Because the stitching process is manual, in some cases, like the 1973 8-inch floppy disc, we need some assistance. We learned how to obtain this assistance from a python script that aligns images and calculates the offset in meters between the stitched images. By combining all the methods together we were able to stitch the MFM data from the floppy disc, which then was added to our on-going effort to measure bit size over the years.
I hope you find this useful, entertaining, and try it yourselves. Please let me know if you use some of this, and as usual, if you have suggestions or requests, don't hesitate to contact me.