[Insight-users] to divide images
Zachary Pincus
zpincus at stanford.edu
Thu, 29 Apr 2004 10:53:18 -0700
--Apple-Mail-7--272700680
Content-Transfer-Encoding: 7bit
Content-Type: text/plain;
charset=US-ASCII;
format=flowed
Hi,
I've modified the itk::ConnectedComponentImageFilter to return a vector
of "submasks" -- that is, images corresponding to each connected blob.
If you want to know where the image belonged in the original, the index
of the region of the returned image will tell you where that sub-mask
came from in the original. This code is attached.
This code is NOT in ITK form -- it's not a filter, or even a class of
any kind. Just a hunk of procedural code. If people want, I could try
to make it into an ImageToImage filter subclass for submission to ITK.
Is this something that would be valuable to others? Let me know.
Zach Pincus
Department of Biochemistry and Program in Biomedical Informatics
Stanford University School of Medicine
PS. I also wrote some code that takes a binary mask and extracts the
region of an image corresponding to the mask (even if the mask is
smaller than the image), optionally rotating it to the principal axes
of the mask. Let me know if you want this too.
--Apple-Mail-7--272700680
Content-Transfer-Encoding: 7bit
Content-Type: text/plain;
x-mac-type=54455854;
x-unix-mode=0644;
x-mac-creator=522A6368;
name="GetConnectedMasks.h"
Content-Disposition: attachment;
filename=GetConnectedMasks.h
/* GetConnectedMasks.h :: Zachary Pincus 3/31/2004
* ----------------------------------------------------
*
*/
#ifndef _GetConnectedMasks_h_
#define _GetConnectedMasks_h_
#include <vector>
namespace CellExplorer { namespace CellExtractor {
template<class TImage> void GetMasks(typename TImage::Pointer maskImage,
std::vector<typename TImage::Pointer > &submasks, bool fullyConnected = false);
}} // end of namespaces
#include "GetConnectedMasks.txx"
#endif
--Apple-Mail-7--272700680
Content-Transfer-Encoding: 7bit
Content-Type: application/text;
x-mac-type=54455854;
x-unix-mode=0644;
x-mac-creator=522A6368;
name="GetConnectedMasks.txx"
Content-Disposition: attachment;
filename=GetConnectedMasks.txx
#ifndef _GetConnectedMasks_txx_
#define _GetConnectedMasks_txx_
#include "GetConnectedMasks.h"
#include <iostream>
#include <map>
#include "itkSmartPointer.h"
#include "itkConnectedComponentImageFilter.h"
#include "itkImageRegionIterator.h"
#include "itkImageRegionConstIterator.h"
#include "itkNumericTraits.h"
#include "itkEquivalencyTable.h"
#include "itkConstShapedNeighborhoodIterator.h"
#include "itkConstantBoundaryCondition.h"
namespace CellExplorer { namespace CellExtractor {
template<class TImage> void GetMasks(typename TImage::Pointer maskImage,
std::vector<typename TImage::Pointer > &submasks, bool fullyConnected) {
typedef TImage ImageType;
typedef typename ImageType::PixelType PixelType;
typedef typename ImageType::SizeType SizeType;
typedef typename ImageType::IndexType IndexType;
typedef typename ImageType::RegionType RegionType;
typedef unsigned short LabelPixelType;
typedef itk::Image<LabelPixelType, 2> LabelImageType;
RegionType masterRegion = maskImage->GetLargestPossibleRegion();
// create an image to hold component labels
LabelImageType::Pointer labelImage = LabelImageType::New();
labelImage->SetRegions(masterRegion);
labelImage->Allocate();
labelImage->FillBuffer(itk::NumericTraits<LabelPixelType>::Zero);
// create an equivalency table
itk::EquivalencyTable::Pointer eqTable = itk::EquivalencyTable::New();
LabelPixelType label, originalLabel, neighborLabel;
LabelPixelType maxLabel = itk::NumericTraits<LabelPixelType>::Zero;
const LabelPixelType maxPossibleLabel = itk::NumericTraits<LabelPixelType>::max();
// Set up the boundary condition to be zero padded (used on label image)
itk::ConstantBoundaryCondition<LabelImageType> BC;
BC.SetConstant(itk::NumericTraits<LabelPixelType>::Zero);
// Neighborhood iterator. Let's use a shaped neighborhood so we can
// restrict the access to face connected neighbors. This iterator
// will be applied to the output image
typedef itk::ConstShapedNeighborhoodIterator<LabelImageType> NeighborhoodIteratorType;
SizeType kernelRadius;
kernelRadius.Fill(1);
NeighborhoodIteratorType nit;
nit = NeighborhoodIteratorType(kernelRadius, labelImage, masterRegion);
nit.OverrideBoundaryCondition(&BC); // assign the boundary condition
// only activate the indices that are "previous" to the current
// pixel and face connected (exclude the center pixel from the
// neighborhood)
//
unsigned int d;
NeighborhoodIteratorType::OffsetType offset;
if (!fullyConnected){
// only activate the "previous" neighbors that are face connected
// to the current pixel. do not include the center pixel
offset.Fill(0);
for (d=0; d < ImageType::ImageDimension; ++d){
offset[d] = -1;
nit.ActivateOffset(offset);
offset[d] = 0;
}
} else {
// activate all "previous" neighbors that are face+edge+vertex
// connected to the current pixel. do not include the center pixel
unsigned int centerIndex = nit.GetCenterNeighborhoodIndex();
for (d=0; d < centerIndex; d++) {
offset = nit.GetOffset(d);
nit.ActivateOffset(offset);
}
}
// along with a neighborhood iterator on the input, use a standard
// iterator on the input and output
itk::ImageRegionConstIterator<ImageType> it;
itk::ImageRegionIterator<LabelImageType> oit;
it = itk::ImageRegionConstIterator<ImageType>(maskImage, masterRegion);
oit = itk::ImageRegionIterator<LabelImageType>(labelImage, masterRegion);
// Mark the output image as either background or unlabeled
it.GoToBegin();
oit.GoToBegin();
while (!it.IsAtEnd()) {
if (it.Get() != itk::NumericTraits<PixelType>::Zero) {
// mark pixel as unlabeled
oit.Set(maxPossibleLabel);
}
++it;
++oit;
}
// iterate over the image, labeling the objects and defining
// equivalence classes. Use the neighborhood iterator to access the
// "previous" neighbor pixels and an output iterator to access the
// current pixel
nit.GoToBegin();
oit.GoToBegin();
while ( !oit.IsAtEnd() ) {
// Get the current pixel label
label = oit.Get();
originalLabel = label;
// If the pixel is not background
if (label != itk::NumericTraits<LabelPixelType>::Zero) {
// loop over the "previous" neighbors to find labels. this loop
// may establish one or more new equivalence classes
typename NeighborhoodIteratorType::ConstIterator sIt;
for (sIt = nit.Begin(); !sIt.IsAtEnd(); ++sIt) {
// get the label of the pixel previous to this one along a
// particular dimension (neighbors activated in neighborhood iterator)
neighborLabel = sIt.Get();
// if the previous pixel has a label, verify equivalence or
// establish a new equivalence
if (neighborLabel != itk::NumericTraits<LabelPixelType>::Zero) {
// if current pixel is unlabeled, copy the label from neighbor
if (label == maxPossibleLabel) {
// copy the label from the previous pixel
label = neighborLabel;
}
// else if current pixel has a label that is not already
// equivalent to the label of the previous pixel, then setup
// a new equivalence. note the use of Lookup() and not
// RecursiveLookup(). this is possible since we keep the
// equivalence table flat.
else if ((label != neighborLabel)
&& (eqTable->Lookup(label) != eqTable->Lookup(neighborLabel))) {
eqTable->Add(label, neighborLabel);
// if we keep the equivalency table up to date, then we
// can use straight calls to Lookup() instead of
// RecursiveLookUp(). This works out to be 3X faster.
eqTable->Flatten();
}
}
}
// if none of the "previous" neighbors were set, then make a new label
if (originalLabel == label) {
// create a new entry label
if (maxLabel == maxPossibleLabel) {
std::cerr << "Number of labels exceeds number of available labels for the output type." <<std::endl;
} else {
++maxLabel;
}
// assign the new label
label = maxLabel;
}
// Finally, set the output pixel to whatever label we have
if (label != originalLabel)
{
oit.Set( label );
}
}
// move the iterators
++nit;
++oit;
}
// Flatten the equavalency table
eqTable->Flatten();
// remap the labels and add the indices of each labeled point to a
// vector that corresponds to each label
typedef std::map<LabelPixelType, std::vector<IndexType> > BlobTableType;
BlobTableType blobTable;
typedef itk::ImageRegionConstIteratorWithIndex<LabelImageType> IndexItType;
IndexItType labelIt(labelImage, masterRegion);
for(labelIt.GoToBegin(); !labelIt.IsAtEnd(); ++labelIt) {
label = labelIt.Get();
// if pixel has a label, write out the final equivalence
if (label != itk::NumericTraits<LabelPixelType>::Zero) {
LabelPixelType flatLabel = eqTable->Lookup(label);
// this is tricky: if gets the vector associated with a given label
// and IF the vector doesn't exist, the [] operator creates a new one!
blobTable[flatLabel].push_back(labelIt.GetIndex());
// Might be able to make this run faster by using an insert method that
// hints where in the map we insert. (always the end!)
}
}
// Now, iterate through the vectors and create a small image from each
for(typename BlobTableType::iterator btIt = blobTable.begin();
btIt != blobTable.end(); ++btIt) {
std::vector<IndexType> indices = (*btIt).second;
// first compute the bounding box of the indices in each vector
// (technically this step could be folded into the act of filling the vectors above)
typedef typename IndexType::IndexValueType IVType;
IVType xMin = itk::NumericTraits<IVType>::max();
IVType yMin = itk::NumericTraits<IVType>::max();
IVType xMax = itk::NumericTraits<IVType>::min();
IVType yMax = itk::NumericTraits<IVType>::min();
for(typename std::vector<IndexType>::const_iterator idxIt = indices.begin();
idxIt != indices.end(); ++idxIt) {
IndexType idx = *idxIt;
if (idx[0] < xMin) xMin = idx[0];
if (idx[0] > xMax) xMax = idx[0];
if (idx[1] < yMin) yMin = idx[1];
if (idx[1] > yMax) yMax = idx[1];
}
// now create an image and put "one" values in every index listed in the vector
SizeType size = {xMax - xMin + 1, yMax - yMin + 1};
IndexType start = {xMin, yMin};
RegionType region;
region.SetSize(size);
region.SetIndex(start);
typename ImageType::Pointer outputImage = ImageType::New();
outputImage->SetRegions(region);
outputImage->Allocate();
outputImage->FillBuffer(itk::NumericTraits<PixelType>::Zero);
for(typename std::vector<IndexType>::const_iterator idxIt = indices.begin();
idxIt != indices.end(); ++idxIt) {
outputImage->SetPixel(*idxIt, itk::NumericTraits<PixelType>::One);
}
// Finally, add that image to the vector that we are passed in.
submasks.push_back(outputImage);
}
}
}} // end namespaces
#endif
--Apple-Mail-7--272700680
Content-Transfer-Encoding: quoted-printable
Content-Type: text/plain;
charset=ISO-8859-1;
delsp=yes;
format=flowed
On Apr 29, 2004, at 10:20 AM, David Llanos wrote:
> Hi all,
> =A0
> I am looking for classes and methods to divide a binary image in =20
> several binary images, in function of the I number of objects that the =
=20
> original binary image has. For example, of the image binar.png, to =20
> obtain four images like uno.jpg.
> I have gotten the first requirement for this objective that is to =20
> label this objects, by means of the following code:
> ...
> I have been looking for some relabeller or labeller method=A0 (in =20
> http://www.itk.org/Insight/Doxygen/html/=20
> classitk_1_1RelabelComponentImageFilter-members.html and in =20
> http://www.itk.org/Insight/Doxygen/html/=20
> classitk_1_1ConnectedComponentImageFilter-members.html) to go keeping =20=
> the different objects of the image in pointers, either to write files =20=
> or to work with them for separate, but I have not found anything with =20=
> that that to be able to work.
> does some form exist of using these classes for my objective, or do I =20=
> have to use another class? In that case, What relabeller or labeller =20=
> method=A0 can I use to relate an object labeled with another class?
> =A0
> Thanks inadvange and regards,
> =A0
> David.
--Apple-Mail-7--272700680--