Interpolate/fill in missing cells in truth-value grid in MATLAB

Question

Consider the 37x101 matrix below:

Each black cell has value 1, the rest of the cells value 0. I would like to fill in the "gaps" by means of cubic spline interpolation, as well as scaling the y axis from 37 to 181. The latter can be done by using the interp1 function, as in:

interp1(37,matrix,181,'pchip')

However, the outcome is now interpolated along the y-axis, but the gaps remain:

I don't want to interpolate along the x-axis, because I want the final matrix to have dimension 181 x 101. I just want to fill in the gaps using the existing cells (of the 181 x 101 matrix).

How can the original matrix (top), be scaled from 37 x 101 to 181 x 101 (without the "smoothing" in the second image), as well as filling in the gaps using some kind of spline interpolation as if this was a proper function?

Answer 1

It appears that your truth value grid has a single one where the true value is in each row. If the true/1 values do in fact create a line through the image, I would recommend parametrize the line with respect to t so that y = fy(t) and x = fx(t). If you're not familiar with this you can find some parametrization info on youtube tutorials or google. The main idea is that if you have , say a truth table that looks like this:

Then you could plot the the location of each pixel with respect to another variable, t and then use interp1(...) on each of these individually. In my case I defined the x and y values as follows:

n = 32;
rand('seed', 1982);
y_orig = 1:n;
x_orig = ceil(n*sin(y_orig/n*pi));

So I can plot as:

t1 = linspace(0,1, n);
plot(t1,x_orig, 'r', 'linewidth', 3);
hold all
plot(t1,y_orig, 'b', 'linewidth', 3);
legend('X', 'Y')

Note that I can get any truth value I want just by using interp1 like this (if you wanted to find the value half way between the 5th and 6th row):

 desiredY = 5.5;
 t= 1:n;
 truthValue= interp1(t, x_orig, desiredY, 'cubic')

But we are looking to make a new image so I chose a more convenient parametrization of t between zero and one. Unfortunately, you may not have x and y off hand, so we need to pull them out of the image. Assuming you have a single true/1 value in each row we can yank out the values with max(...):

[maxVals, x1] = max(data,[],2);
x1(maxVals == 0) = [];
y1 = find(maxVals ~= 0);

Some form of find on each row would also work. If you have a truth value in each row then y1 should equal 1:n. The max function returns the index of the max in dimension 2 in the second return value. I use the next two lines to remove any entries where there truth table was empty (max is zero) and then y1 = 1:n minus those entries that were empty.

A quick and dirty way to get lots of points along this line is:

t2 = linspace(0,1,1024);
x2 = interp1(t1, x1, t2, 'cubic');
y2 = interp1(t1, y1, t2, 'cubic');

I can then plot the original points/image and this newly discovered finer line together like this:

imagesc(data);
hold all;
plot(x2,y2, 'linewidth', 2);
axis image
colormap(flipud(colormap(gray)));

To get this:

Finally, you can quickly turn this into a new image by scaling the parametrization up. My method is not particularly efficient for clarity:

y2_scaled = floor((y2(:)-1)*scaleValue + 1);
x2_scaled = floor((x2(:)-1)*scaleValue + 1);

scaleValue = 2;
data2 = zeros(n*scaleValue);
for ind = 1:length(x2_scaled)
    data2(y2_scaled(ind),x2_scaled(ind)) = 1;
end

Which results in:

Note that this table has connected all the points and you now have multiple true/1's in each row. This is because I chose a very small step size for t2. You could fix this by either choosing t2 smarter, skipping multiple values in each row, or average the location of each indices in each row. Or ignoring this issue.

To fix t2 with the scaling value, you could use

t2 = linspace(0,1,n*scaleValue);

to get only one true/1 per row in the above code.

Also, if you want to only scale one dimension, you could do it like this:

y2_scaled = floor((y2(:)-1)*scaleValue + 1);
x2_scaled = floor((x2(:)-1) + 1);

scaleValue = 2;
data2 = zeros(n*scaleValue,n);
for ind = 1:length(x2_scaled)
    data2(y2_scaled(ind),x2_scaled(ind)) = 1;
end

Answer 2

I see this as a bitmap, so why not a clamped blur?

I=yourmatrixhere

%gaussian blur
%    it looks like bump-to-bump distance is 3 empties
%    therefore hsize should be about 7
%    it looks like bump vertical size is about 4
%    therefore simga should be about 10

hsize=[3 3];
sigma = 10;
h=fspecial('gaussian',hsize,sigma)

I2=imfilter(I,h,'replicate');

At this point you have spread information to adjacent columns, but you need to "tidy up" from continuous to binary.

%threshold
th = 0.25;

I3=zeros(size(I));
ind=find(I>=th);
I3(ind)=1;

At this point, I3 is your matrix of interest to do the "erode" or interpolation.