In part 1, we have discussed the steps to find the tiling from the scanned paper. Just to remind you, there were:

- Rotate the image to be orthogonal to x-y axis
- Find the grid
- For each edge on the grid, decide whether it is “black edge” or “no edge”
- Break the grid into several parts, and process each one individually.
- Somehow, by using the prior knowledge on the 12 Pentomino pieces, find where they are

One small thing, I did show the code of detecting the angle, but did not write the actual rotation command. It is as simple as:

rotatedIm = imrotate(clippedIm,90-angle);

I am skipping a small part of my code that removes white background. We can go back to it later. Here are two images, before clipping and rotation, and after:

You can see that now the lines of the grid are truly vertical.

Now, let’s focus on the second task – “Finding the grid”. Just to be clear, finding the grid means finding a vector of X and a vector of Y that represent the places where the lines in the image are. The most obvious choice for detecting periodic patterns in image processing is DFT (Discrete Fourier Transform), which is called FFT in Matlab. I am talking about summing up the image in either x or y direction, and using `fft`, not using `fft2`, which is not needed after we found the angle and rotated.

FFT can only be used as initial guess, since it is not exactly accurate. The reason is due to its quantization. For instance, imagine that the distance between the grid lines is about 20 pixels. The size of the image in x direction is 1024. Can you detect it? Well, FFT checks the signal decomposition into several frequencies. We need to find the frequency of 1024/20 = 51.2 , which is between 51 and 52. So how accurate is our estimation? We will either find the maximum correlation between the signal and the corresponding sine at the frequency of 51 or 52. Simple math shows that 1024/51 = 20.0784, 1024/51 = 19.6923. That means that we will either have a mistake of 0.784 or 0.307 pixels. That doesn’t sound much, right? After all, something less than a pixel must be small. Wrong! Don’t forget that a small mistake in the the distance is multiplied by the amount of lines, that is about 50. That can easily become as much as 0.784*50 = 39.20 pixels. Now that is a **huge error.**

Let’s see what happens just for the fun.

function cols = DetectColumnsWithFFT(im) s = sum(im,1); f = fft(s); magnitude = abs(f); phase = angle(f); %Remove the DC magnitude(1) = 0; [~,maxFreqIndex] = max(magnitude); %Estimate distance between grid elements gridLength = size(im,2)/ maxFreqIndex; gridPhase = phase(maxFreqIndex); startIndex = 1 + (gridLength / 2) * (gridPhase/pi); cols = startIndex: gridLength: size(im,2); end

Here is the result:

The information from the phase allowed us to find the first column exactly. Yet the frequency cannot be found precisely, and eventually the error starts to accumulate. To sum up, FFT is valuable because it is very fast, robust and gives a good *initial* estimate. The main problem is that by using it, we cannot find the frequency accurately enough

The way to solve this issue is to look for each of the grid lines independently. If I sum the image row-wise, I will get something like this:

The valleys correspond to the grid lines. I am going to use an ad-hoc heuristical solution. First, find a middle line and threshold by its value

Then, I get a square-like wave that consists of [0,1]. I can find the center of each of the zeros sequence. This will give me almost all of the lines.

There are some missing lines! If you take a look again at the graph above, you will understand why. I will show you next time how to fill them.

Here is the code.It is the same both for rows and columns detection. If I want to find columns, the `dim` parameter should be 1, and in case of rows 2. Can you find out what that small trick that I am using in lines 3-5 does? I will leave it as as something to think about.

function centers = FindFrequencyOnStraight(greenIm,dim) peakWave = sum(greenIm,dim); places = linspace(1,numel(peakWave),10); places = round(places); peakWaveClipped = peakWave( places(2):places(end-1)); nMax = max(peakWaveClipped); nMin = min(peakWaveClipped); nMiddle = (nMax + nMin)/2; bThresholdedPeakWave = peakWave>nMiddle; centers = regionprops(~bThresholdedPeakWave,'Centroid'); centers = [centers.Centroid]; if dim==1 centers(2:2:end) = []; else centers(1:2:end) = []; end end

Next time I am going to show how to find those missing lines, and explain about how to find out the edges of the pieces.

**Once again, I would like to remind you that this is one of my first posts and a feedback is very welcome!**

Way to go! Found it very interesting, even that this is totally not my area of expertise.