[Insight-users] to divide images

David Llanos gva02@elai.upm.es
Mon May 10 18:29:22 EDT 2004


This is a multi-part message in MIME format.

------=_NextPart_000_0042_01C436C5.188F2CB0
Content-Type: text/plain;
	charset="iso-8859-1"
Content-Transfer-Encoding: 7bit

Hi Zach,

I have been trying to use your method GetMasks in all the ways, but I am not
able to get it.
Please, I request you that you help me to correct the code that I attach you
next, or that you send me an example of as using your bookstore and the
method:
----------------------------------------------------------------------------
----
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkImage.h"
#include "itkConnectedComponentImageFilter.h"
#include "itkRelabelComponentImageFilter.h"
#include "itkImageMomentsCalculator.h"

#include "GetConnectedMasks.h"

int main ( int argc, char * argv[] )
{
 typedef unsigned char   PixelType;
 const unsigned int   dimension = 2;
 typedef itk::Image<PixelType,dimension>  ImageType;
 typedef itk::ImageFileReader<ImageType>   itkReaderType;
 typedef itk::ImageFileWriter<ImageType>   itkWriterType;
 typedef itk::ConnectedComponentImageFilter<ImageType,ImageType>
  itkConnectedComponentFilterType;
 typedef itk::RelabelComponentImageFilter<ImageType,ImageType>
  itkRelabelComponentFilterType;

 typedef std::vector<ImageType> &objetos;

 itkReaderType::Pointer reader  = itkReaderType::New();
 reader->SetFileName( argv[1] );
 reader->Update();

 CellExplorer::CellExtractor::GetMasks(reader, &objetos, TRUE);

 std::vector<unsigned long>::const_iterator it;
 int i;
 for(i = 0, it =  objetos.begin(); it != objetos.end(); ++i, ++it) {
  itkWriterType::Pointer submascara = itkWriterType::New();
  submascara->SetInput( submascara->GetOutput() );
  submascara->SetFileName( "submascara.png" );
  submascara->Write();
 }
 return 0;
}
----------------------------------------------------------------------------
-------------------
E:\label\etiquetar.cxx(95) : error C2275: 'objetos' : illegal use of this
type as an expression
        E:\label\etiquetar.cxx(36) : see declaration of 'objetos'
E:\label\etiquetar.cxx(100) : error C2228: left of '.begin' must have
class/struct/union type
E:\label\etiquetar.cxx(100) : error C2228: left of '.end' must have
class/struct/union type
E:\label\etiquetar.cxx(103) : error C2661: 'GetOutput' : no overloaded
function takes 0 parameters
Error executing cl.exe.

etiquetar.exe - 4 error(s), 0 warning(s)
----------------------------------------------------------------------------
---------------------

Also, I attach you the image that I am trying to divide and the bookstores
that you sent me, for if somebody also wants to take them a look. I am
learning studying them a lot.


----- Original Message ----- 
From: "Zachary Pincus" <zpincus@stanford.edu>
To: "David Llanos" <gva02@elai.upm.es>
Sent: Wednesday, May 05, 2004 8:05 PM
Subject: Re: [Insight-users] to divide images


> Here again is the source for my simple masked region extractor.
>
> The header has just one function: GetMasks. You call this function,
> passing in:
> (1) A mask image where the foreground is value 1 and the background is
> value 0
> (2) A vector that the "submasks" will be placed in. Each submask is a
> new image of the same type. The regions set on the submask image tell
> you where those images came from in the original.
> (3) Optionally, a boolean value that tells the function whether to look
> for "face" or "fully" connected blobs. (face = 4-connected, fully =
> 8-connected, in 2D.) If you are looking for wide blobs that cover some
> area, the default of "false" is fine. If you are trying to extract
> single-pixel wide lines, you should use fully connected mode. (That is,
> send "true").
>
> Zach

------=_NextPart_000_0042_01C436C5.188F2CB0
Content-Type: image/png;
	name="binar.png"
Content-Transfer-Encoding: base64
Content-Disposition: attachment;
	filename="binar.png"

iVBORw0KGgoAAAANSUhEUgAAAecAAAFsCAAAAADTES6TAAAACXBIWXMAAC4jAAAuIwF4pT92AAAE
bUlEQVR4nO3dy1ZVMRBFUf7/p7EhooLIeSSVnVNzdmyyU0s6cMfg5QUAAAAAAAAAAAAAAAAAAAAA
AAAAAAAAAAAAAADo6vWD1XuY42NnpR/pc2axH+Utp86P9lVenZ/k28o6P8GBzELv71Bmnbd3rLPU
uzvcWemdncgs9K5ORdZ6P1cC67wfnXu40VnojdzprPQ+dO5B5x507uFeZ6l3ofPj3U6s8xZ07kHn
HnTuYUxnocMNyqx0uHGdpU42tLPQsXTuQecWxmbWOZXOLcjcgsw96NyDzj3o3IPOPajchM496NyD
zj2M7Cx0Lp170LkHnXvQuQeZe9C5CZ170LkHnXvQuQuZm9C5B5170LkHmXuQuAedu5C5B5270LkH
nXuQuQed29C5EZF70LkHnRtRuQeRe9C5B5170LiL954f8sr8ML9yfvw2lvmhfAMDAAAAAFT69ENZ
P6V9pN9Z/db0yf79azWhn+bLj0HI/Shff+JF50f46+MtOj/WN4V1foBDjXXe3uHOQu/reGSl93Wy
stinBBzrLdi1zmp/J+Jif3zhO52F/lrE3W7F1fmIiMvpPF3E5XSeLuJ4Ok8XcT2d50o535aV379o
/H+znAPuGHr5gKNyDjgmc9Gd1y84J+iEN9LWnnn5gAuCLnixavmZY4acEXTBc1PWnTdmyAlJNzy9
Zc15Y4Ycl3XDC2sWnDdoylFhN7wyp3plzpIToi44tPOsmTFDTsm64QadY4acE3XDoZnnbMxaM313
zpLSjRv/DZmY1TdOWDVx58+sxay+c8MdJo6fUzE9ZEbhxK3/NGfK6nsnrNmo88IdhRs7dk7ZUboy
a03J9JQdhSuz1hSNjxlStzJrTdX4mCFlQ6PGlI2PGVI2NGpM2fiYIVVLk7YUro8ZUrQ0aUvl/Jgh
RUuTtlTOz1lSsjRpS+n+nCUlS5O2VO7PWVIzM25Q1f6cJRUr4waVPSFkRtHKuEFlb4gYUbYyb1HV
GyJGlK3MW1T1hogRZSvzFlW9IWJE2cq8RVVviBhRNTJwUs071i8oXZm3qOIhAROKp6btqXlIwITi
qXGDKh4SMOHIypFbR+5aKmDpwFu+vup89yEBE47tHLt24LDFVo8cd8kZS2OHzXnK4i+/7pyxwy5Z
OXLcJWcsjR121bKR4y45Y2nssKuWTRx3yRlTc5ddtGzduEvO2Bs77IY164Zdcsrg2GE31Y97+0qh
5xw0K67zMqHHHDJryrI9pR5zyK4507aUeswhu+ZM21LqLUfsmrVtR6m3HLFr1rYdpZ5xSGSV36Ue
UubxEs+o8xx5Z9R5hrwzyjxF3B01nijolvc6C/1fhy5YdUmdp4m6oM4TJV1Q5xLLL6hzibd7/fxH
56d7P1v5AXXuQeceVO5B5iZEbkLnRnTuQecedO5C4DZk7kHnHnRu4j+BdX4W38OtaNyDygAAAAAA
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA0NAP9SymMqfm
gH8AAAAASUVORK5CYII=

------=_NextPart_000_0042_01C436C5.188F2CB0
Content-Type: application/octet-stream;
	name="GetConnectedMasks.h"
Content-Transfer-Encoding: quoted-printable
Content-Disposition: attachment;
	filename="GetConnectedMasks.h"

/* GetConnectedMasks.h :: Zachary Pincus 3/31/2004
 * ----------------------------------------------------
 *
 */
=20
#ifndef _GetConnectedMasks_h_
#define _GetConnectedMasks_h_

#include <vector>

namespace CellExplorer { namespace CellExtractor {
//namespace itk
=09
	template<class TImage> void GetMasks(typename TImage::Pointer =
maskImage,
		std::vector<typename TImage::Pointer > &submasks, bool fullyConnected =
=3D false);

}} // end of namespaces

#include "GetConnectedMasks.txx"
#endif
------=_NextPart_000_0042_01C436C5.188F2CB0
Content-Type: application/octet-stream;
	name="GetConnectedMasks.txx"
Content-Transfer-Encoding: quoted-printable
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 {
//namespace itk

template<class TImage> void GetMasks(typename TImage::Pointer maskImage,
		std::vector<typename TImage::Pointer > &submasks, bool fullyConnected) =
{
=09
	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;
=09
	RegionType masterRegion =3D maskImage->GetLargestPossibleRegion();
=09
	// create an image to hold component labels
=09
	LabelImageType::Pointer labelImage =3D LabelImageType::New();
	labelImage->SetRegions(masterRegion);
	labelImage->Allocate();
	labelImage->FillBuffer(itk::NumericTraits<LabelPixelType>::Zero);
=09
	// create an equivalency table
	itk::EquivalencyTable::Pointer eqTable =3D =
itk::EquivalencyTable::New();

	LabelPixelType label, originalLabel, neighborLabel;
	LabelPixelType maxLabel =3D itk::NumericTraits<LabelPixelType>::Zero;
	const LabelPixelType maxPossibleLabel =3D =
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	 =3D 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=3D0; d < ImageType::ImageDimension; ++d){
			offset[d] =3D -1;
			nit.ActivateOffset(offset);
			offset[d] =3D 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 =3D nit.GetCenterNeighborhoodIndex();
		for (d=3D0; d < centerIndex; d++) {
			offset =3D 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 =3D itk::ImageRegionConstIterator<ImageType>(maskImage, =
masterRegion);
	oit =3D itk::ImageRegionIterator<LabelImageType>(labelImage, =
masterRegion);
=09
	// Mark the output image as either background or unlabeled
	it.GoToBegin();
	oit.GoToBegin();
	while (!it.IsAtEnd()) {
		if (it.Get() !=3D itk::NumericTraits<PixelType>::Zero) {
			// mark pixel as unlabeled
			oit.Set(maxPossibleLabel);
		}
	=09
		++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 =3D oit.Get();
		originalLabel =3D label;
	=09
		// If the pixel is not background
		if (label !=3D 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 =3D 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 =3D sIt.Get();

				// if the previous pixel has a label, verify equivalence or
				// establish a new equivalence
				if (neighborLabel !=3D itk::NumericTraits<LabelPixelType>::Zero) {
					// if current pixel is unlabeled, copy the label from neighbor
					if (label =3D=3D maxPossibleLabel) {
						// copy the label from the previous pixel
						label =3D 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 !=3D neighborLabel)
								&& (eqTable->Lookup(label) !=3D 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 =3D=3D label) {
				// create a new entry label
				if (maxLabel =3D=3D maxPossibleLabel) {
					std::cerr << "Number of labels exceeds number of available labels =
for the output type." <<std::endl;
				} else {
					++maxLabel;
				}

				// assign the new label
				label =3D maxLabel;
			}

			// Finally, set the output pixel to whatever label we have
			if (label !=3D 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
=09
=09
	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 =3D labelIt.Get();
		// if pixel has a label, write out the final equivalence
		if (label !=3D itk::NumericTraits<LabelPixelType>::Zero) {
			LabelPixelType flatLabel =3D eqTable->Lookup(label);
		=09
			// 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!)
		}
	}
=09
	// Now, iterate through the vectors and create a small image from each
	for(typename BlobTableType::iterator btIt =3D blobTable.begin();=20
		btIt !=3D blobTable.end(); ++btIt) {
=09
		std::vector<IndexType> indices =3D (*btIt).second;
	=09
		// 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)	=09
	=09
		typedef typename IndexType::IndexValueType IVType;
	=09
		IVType xMin =3D itk::NumericTraits<IVType>::max();
		IVType yMin =3D itk::NumericTraits<IVType>::max();
		IVType xMax =3D itk::NumericTraits<IVType>::min();
		IVType yMax =3D itk::NumericTraits<IVType>::min();
	=09
		for(typename std::vector<IndexType>::const_iterator idxIt =3D =
indices.begin();=20
			idxIt !=3D indices.end(); ++idxIt) {
				IndexType idx =3D *idxIt;
				if (idx[0] < xMin) xMin =3D idx[0];
				if (idx[0] > xMax) xMax =3D idx[0];
				if (idx[1] < yMin) yMin =3D idx[1];
				if (idx[1] > yMax) yMax =3D idx[1];
		}
	=09
		// now create an image and put "one" values in every index listed in =
the vector
		SizeType size =3D {xMax - xMin + 1, yMax - yMin + 1};
		IndexType start =3D {xMin, yMin};
		RegionType region;
		region.SetSize(size);
		region.SetIndex(start);
	=09
		typename ImageType::Pointer outputImage =3D ImageType::New();
		outputImage->SetRegions(region);
		outputImage->Allocate();
		outputImage->FillBuffer(itk::NumericTraits<PixelType>::Zero);

		for(typename std::vector<IndexType>::const_iterator idxIt =3D =
indices.begin();=20
			idxIt !=3D indices.end(); ++idxIt) {
				outputImage->SetPixel(*idxIt, itk::NumericTraits<PixelType>::One);
		}
	=09
		// Finally, add that image to the vector that we are passed in.
		submasks.push_back(outputImage);
	}
}


}} // end namespaces



#endif
------=_NextPart_000_0042_01C436C5.188F2CB0--




More information about the Insight-users mailing list