Automatically trim/remove/crop black borders/margins from images / volumes

Is there any way to easily and efficiently automate "trimming / removing / cropping" of "black / 0", "borders / margins" from "images / volumes" in Matlab?

In particular, I am working with 3D volumes and even 4D volumes of medical image datasets, in which the CT or MR reconstruction algorithms occasionally leave a geometric, contiguous, unused, black margin around the images.

For purposes of doing image analysis and statistics, it helps to remove this useless, spurious, zero data which tends to be problematic and contaminating when performing image analysis.

[Eric Diaz]

Is anybody else available to help with this problem?

[Benjamin Avants]

Is the border region always the same size for a given type of data?

[Eric Diaz]

I am unable to say for certain whether the black, border region is always the same, as it is an artifact of the reconstruction process. Thus, I think it would be safer to make a more generalizable algorithm.

[Eric Diaz]

Can anyone else help with this problem? It seems that everyone yet has given up or is no longer interested in helping.

[Image Analyst]

Or is just too busy at the moment to spend the time needed.

Answer 1

[Image Analyst]

Get the black pixels imageArray == 0. Then invert and call imclearborder and subtract or XOR the two to get only those black pixels touching the border. Then use that mask to set the pixels to the desired gray level imageArray(mask) = desiredGrayLevel.

[Eric Diaz]

Thanks for the suggestion. I will give that a try and let you know if I have any problems.

[Eric Diaz]

well... the first curiosity I've run into is that after running the following commands

blackPixels = img==0;
whitePixels = ~blackPixels;
im2         = imclearborder(whitePixels);
mask        = xor(blackPixels,im2);

 size(img)

ans =

   320   320    10     5

 size(im2)

ans =

   320   320    10     5

The size of my image array ('img') and "imClearBorder" array ('im2') are identical, which was not what I was expecting. Is this supposed to work that way?

[Eric Diaz]

Oh, wait, I'm not sure what you mean by set the img(mask) to desired gray level.

Would this actually remove those pixels and reduce the matrix size or just set them to a particular value.

If it is the latter, I'm not sure I understand the value of doing that since I already have a known value for those pixels, namely 0.

[Image Analyst]

I don't have any data to work with. The desired gray level is the gray level that you want the outer black pixels to be instead of 0. Maybe it's 20 or 100 or 222 - I have no idea what you want to do with them.

[Eric Diaz]

Ah, okay, that's what I was thinking you meant.

Well, I was actually preferring to just remove them from the matrix, rather than substitute a value in for the zero.

Any ideas?

[Image Analyst]

You cannot. Arrays must be rectangular. Period. They cannot have ragged edges.

[Eric Diaz]

Right. I knew that. I was trying to think of a way to create the largest possible array without any of the zeros. So, in other words, some real data will have to be thrown out. However, I would prefer to throw out some real data near the edges of the image, which is usually background, rather than retain useless zeros or substitute an arbitrary value into those zeros. Both of the latter options tend to be problematic for analysis, whereas the former does lose some data but in the minimized sense (i.e. - largest possible array after exclusion of junk pixels).

[Image Analyst]

You can crop if you want but that's usually not needed. More common is masking to just tell it what pixels to ignore.

[Eric Diaz]

Ah, yes, good thinking. Not sure why that hadn't occurred to me.

I just have a quick question though, which I fear may be the reason I am hesitant to try that approach.

Let's say I have this 3D volume and I perform the following,

blackPixels = img==0;
whitePixels = ~blackPixels;

I now have a can use the whitePixels as my include mask and the blackPixels as my exclude mask.

If my next steps include the following actions:

 1) reshape my 3D image volume into a vector (ImgVec)

 2) perform an unsupervised K-Means clustering, returning a vector (clusterIndx) containing indices labeled by cluster.

 3) reshape clusterIndx to 3D volume with labels.

I am trying to figure out how I might apply the mask that without screwing up the reshaping transform from 3D --> 1D --> 3D. Any thoughts?

For example, one of these 1D vector would be fed into the clustering algorithm. But then I would need to retransform that to 3D. However, due to differences in size (length), the transform to 3D cannot be made.

maskedImg = img(whitePixels);

imgVec=img(:);

size(maskedImg)

ans =

     4991735           1

 size(imgVec)

ans =

     5120000           1

[Eric Diaz]

Can anyone help with this?

[Jurgen]

So you selected your pixels of interest with a mask, now you want to give them nominal values (e.g. integers 1-10) based on your clustering result. Your problem is that you don't know how to map the vector back to the pixels of interest?

You already have the mask...so replace the 1 values in the mask with the labels to get a label mask. Then apply the new label mask to img.

[Eric Diaz]

I think you understand my problem correctly.

I think I understand your explanation but am not sure. I will give it a try.

Thank you for your suggestion.

[Image Analyst]

What are you sending into kmeans()? If you turn your 3D binary volume into a column major vector, then all you'll have is 0's and 1's. What kind of information do you expect to be in clusters? What do you want to do anyway? Do you want to label voxels into clusters, even though they may not be touching? If so, you're going to have to do kmeans on the (x,y,z) coordinates, not the true/false (1/0) value of the voxels.

[Eric Diaz]

Well, I was going to send in the masked image vector

 maskedImgVec = img(whitePixels);

in place of what I used to send in, which was just the image vector

 imgVec = reshape(img,x*y*z,1); % 3D --> 1D

I presume, based on your comment, that doing the reshape command into a column major vector is different from the 'colon index call', i.e., img(:)?

Based on prior work I have done, sending the reshaped 3D volume column major vector into kMeans has produced results consistent with what I was wanting.

The results were voxels labeled into clusters in a column major format.

 [cluster_index, cluster_cntr] = kmeans(imgVec,numK,'distance','sqEuclidean','Replicates',5,'Options',opts)

Those cluster_index results were then reshaped back into 3D.

voxel_labels = reshape(cluster_index,x,y,z); % 1D --> 3D

When I would show a slice from the 3D labeled voxels

figure;imshow(voxel_labels(:,:,5),[]);

I would get a nice result consistent with what I was expecting. Labeled areas of my image.

My problem right now is that in the initial 3D volume there are often useless, "junk" pixels that are an artifact of reconstruction. When reading in the 12 bit images into a UINT16 matrix, there is a sometimes a geometric, congruent rim of pixels around the border which =0.

If I feed these pixels into kMeans, it will also try to cluster these as another cluster, diminishing the degrees of freedom it has to cluster the real data.

However, my other problem now is that if I use a masked image vector,

maskedImgVec = img(whitePixels); % 3D --> 1D

and feed that 1D maskedImgVec into kMeans, I don't know if the results, cluster_index, can be remapped 1D --> 3D, into 3D voxel_labels, because of the different lengths of the vectors, as I showed earlier.

[Eric Diaz]

@Jurgen I was not able to get that technique to work. Perhaps I misunderstood your suggestion. Would you mind re-explaining your suggestion with more detail? Thank you.

[Image Analyst]

This is getting too long, involved, and complicated. Can you upload a very small 3D volume with a small snippet of code so we have something not-so-abstract to work with?

[Eric Diaz]

Sure...I tried my best to explain, but I guess sometimes just playing around with some data works best. Obviously, none of things I am still trying to find out have changed.

load('img.mat')
[rows cols nzs] = size(img);
figure;imshow(img(:,:,5),[])

blackPixels = img==0;
whitePixels = ~blackPixels;
size(blackPixels);

figure;imshow(blackPixels(:,:,5),[])
figure;imshow(whitePixels(:,:,5),[])

ImgVec = reshape(img,rows*cols*nzs,1);
size(ImgVec);
maskedImgVec = img(whitePixels);
size(maskedImgVec);

opts =  statset('Display','iter');
[cluster_idx1, cluster_center1] = kmeans(ImgVec,5,'distance','sqEuclidean','Replicates',5,'Options',opts);

[cluster_idx2, cluster_center2] = kmeans(maskedImgVec,5,'distance','sqEuclidean','Replicates',5,'Options',opts);

voxel_labels = reshape(cluster_idx1,rows,cols,nzs);
figure;imshow(voxel_labels(:,:,5),[])

[Benjamin Avants]

Wow, those borders are troublesome.

I can think of a few things that might help.

One option is to simply set all the 0 valued pixels to NaNs and use functions that disregard NaNs for analysis.

Another option would work if the valuable data in the images is never closer to the edge of the image than the thickest part of the border. You could either find (for each new image) or hard code the number of pixels in the thickest part of the border and simply extract the center of the image.

Let's say the thickest part of the border is 50 pixels:

newImg = img(50:end-50,50:end-50,:);

Matlab can only work with regularly shaped arrays so you will either have to crop each edge by the thickest border on that edge or convert the values of the border pixels to a value that MATLAB will ignore or that won't interfere with your image processing.

[Eric Diaz]

Thanks Ben, for your input.

Yeah, I agree that the borders are challenging. It's kind of what makes it an interesting problem.

Ideally, the solution would completely preserve the "valuable" data while eliminating the "junk" data.

I have tried the hard coding solution several times when playing around with the data, just as you wrote it ... and while it works, it's not the same as finding an "elegant" solution, if you know what I mean.

[Benjamin Avants]

Completely. I hate using work arounds and unoptimized solutions unless I really have to.

[Image Analyst]

Eric: See if the attached demo, test.m below the image, does what you want. Sorry for the delay but I was swamped with work at my real job.

[Eric Diaz]

Thank you Image Analyst for your help.

I am sorry if I was impatient. I thought you had abandoned me, when I didn't hear back after the email.

Your solution appears very elegant and near ideal.

My only issue with the solution is that it actually trims too much data. The reason that this is important is because that while the purely "zero value" data around the rim is actually "junk", the close to zero, but not actually zero, dark background data is "not junk" and important to preserve. The non-zero background data is used to help estimate the noise in the image.

[Image Analyst]

Then change the thresholding line so that you get everything except pure zero:

binaryImage = img > 0;

[Eric Diaz]

Very Elegant & Ideal solution! Thank you very much Image Analyst!

[Image Analyst]

You haven't officially "Accepted" yet, so what's left that you need my help on? If nothing then go ahead and Accept the answer please.

[Eric Diaz]

Thanks for the reminder!

[Eric Diaz]

I just wanted to add that after further investigation, the solution turned out to be not quite "ideal" but it was still a good solution. The solution did not quite completely eliminate "junk data" while maximizing retained data. Rather, it minimized junk data and maximized retained data.

The solution was still very helpful as it did help give me a springboard of ideas to hack upon and try to perfect. I think I am close to tweaking the solution to get my ideal.

[Image Analyst]

I thought my original solution wasn't bad. I don't know why you want pixels with values less than 300 anyway. What is junk? Can you define it based solely on gray level? Solely on location, like outside the main body? Why do you want to crop anyway? What do you want to find, assuming you could do what you want and get rid of junk data?

[Eric Diaz]

The "junk" is the data that is 0. It is an artifact in more ways than one.

The data that is not zero, but with very small grayscale value is actually useful. There are usually thousands of pixels with very small grayscale values which can be used to estimate the noise.

If noise estimation wasn't important, there wouldn't be a need to preserve anything below a certain threshold, but that's not the case.

The junk is always outside the main body. Unfortunately, there are sometimes digitization artifacts within the body, i.e. - sparsely scattered zero pixels. As these pixels are within the main body, they must be retained as pertinent digital artifacts.

The idea is that the noise estimation needs to be as robust and accurate as possible, because when trying to estimate the true signal, an accurate estimate of the noise is needed.

[Image Analyst]

So you want to get all low value pixels outside the main body that are not exactly zero? If so, that's really easy. Just threshold at something like 300 and at zero and OR them and invert then extract

mask = grayImage > 300; % Mask is bright stuff.
% Fill in the body
mask = imfill(mask, 'holes'); % Mask is whole solid body.
% OR it in with the zeros
mask = mask | (grayImage == 0); % Mask now includes pure zeros.
% Extract pixels that are not masked
darkNonZeroOutsidePixels = grayImage(~mask);

[Eric Diaz]

That is a solution that works. I like it.

I think I have just figured out how to fuse the two solutions to beta solution which works best for me.

I can use the aggressive crop to retain main body and perform means clustering on that subset of images. And I can use the above, XOR solution to perform noise estimation on just the useful pixels.

The additional problem, which I needed to address was the degrees of freedom when performing clustering.

Thanks for your help and suggestions, Image Analyst.

[lt c]

mask = grayImage > 300; % Mask is bright stuff.

% Fill in the body

mask = imfill(mask, 'holes'); % Mask is whole solid body.

% OR it in with the zeros

mask = mask | (grayImage == 0); % Mask now includes pure zeros.

% Extract pixels that are not masked

darkNonZeroOutsidePixels = grayImage(~mask);

Quick question:

I ran code above and get a column vector (darkNonZeroOutsidePixels). Does it mean to be a part of 2D image (without the middle bright and the 0 value edge)?