KWStyle - itkBloxCoreAtomImage.txx

Matrix View

Description

---

1		/*=========================================================================
2
3		Program:   Insight Segmentation & Registration Toolkit
4		Module:    $RCSfile: itkBloxCoreAtomImage.txx.html,v $
5		Language:  C++
6		Date:      $Date: 2006/01/17 19:15:33 $
7		Version:   $Revision: 1.4 $
8
9		Copyright (c) Insight Software Consortium. All rights reserved.
10		See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
11
12		This software is distributed WITHOUT ANY WARRANTY; without even
13		the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
14		PURPOSE.  See the above copyright notices for more information.
15
16		=========================================================================*/
17		#ifndef __itkBloxCoreAtomImage_txx
18		#define __itkBloxCoreAtomImage_txx
19
20		#include "itkBloxCoreAtomImage.h"
21
22		#include "itkImageRegionIterator.h"
23		#include "itkImageRegionConstIterator.h"
24		#include "itkConicShellInteriorExteriorSpatialFunction.h"
25		#include "itkEllipsoidInteriorExteriorSpatialFunction.h"
26		#include "itkFloodFilledSpatialFunctionConditionalIterator.h"
27
28		#include "vnl/vnl_matrix.h"
29
30		namespace itk
31		{
32
33		template <unsigned int dim>
34		BloxCoreAtomImage<dim>
35		::BloxCoreAtomImage()
36		{
37		m_MedialNodeCount = 0;
38		m_NodePointerList = new std::vector<BloxCoreAtomPixel<NDimensions>*>();
39		}
40
41		template <unsigned int dim>
42		BloxCoreAtomImage<dim>
43		::~BloxCoreAtomImage()
44		{
45		delete m_NodePointerList;
46		}
47
48		template <unsigned int dim>
49		void
50		BloxCoreAtomImage<dim>
51		::DoEigenanalysis()
52		{
53		itk::ImageRegionIterator<Self> bloxIt =
54	IND	****itk::ImageRegionIterator<Self>(this, this->GetLargestPossibleRegion() );
55
56		for(bloxIt.GoToBegin(); !bloxIt.IsAtEnd(); ++bloxIt)
57		{
58		( &bloxIt.Value() )->DoCoreAtomEigenanalysis();
59		}
60		}
61
62		template <unsigned int dim>
63		void
64		BloxCoreAtomImage<dim>
65		::DoCoreAtomVoting()
66		{
67		//cerr << "BloxCoreAtomImage::DoCoreAtomVoting()" << endl;
68
69		// Iterator to access all pixels in the image
70		ImageRegionIterator<Self> bloxIt =
71	IND	****ImageRegionIterator<Self>(this, this->GetLargestPossibleRegion() );
72
73		// Pointer for accessing pixels
74		BloxCoreAtomPixel<NDimensions>* pPixel = 0;
75
76		// Results of eigenanalysis from each pixel
77		typename BloxCoreAtomPixel<NDimensions>::EigenvalueType eigenvalues;
78		typename BloxCoreAtomPixel<NDimensions>::EigenvectorType eigenvectors;
79
80		// Results of eigenanalysis from each pixel
81		typename BloxCoreAtomPixel<NDimensions>::EigenvalueType sf_eigenvalues;
82		typename BloxCoreAtomPixel<NDimensions>::EigenvectorType sf_eigenvectors;
83
84		unsigned int voterCount = 0;
85		for(bloxIt.GoToBegin(); !bloxIt.IsAtEnd(); ++bloxIt)
86		{
87		// Get a pointer to the pixel
88		pPixel = &bloxIt.Value();
89
90		// If there are no core atoms in this pixel, it doesn't get to vote
91		if( pPixel->empty() )
92		{
93		continue;
94		}
95		//populate the NodePointerList
96		m_NodePointerList->push_back(pPixel);
97
98		voterCount++;
99
100		// Get eigenanalysis results
101		eigenvalues = pPixel->GetEigenvalues();
102		eigenvectors = pPixel->GetEigenvectors();
103
104	LEN	//std::cerr << "eigen values: " << eigenvalues[0] << " " << eigenvalues[1] << " " << eigenvalues[2] << std::endl;
105
106		// Ellipsoid axis length array
107		Point<double, NDimensions> axisLengthArray;
108
109		// Compute first length
110		axisLengthArray[0] = 0.5 * pPixel->GetMeanCoreAtomDiameter();
111
112	LEN	// printf("Mean core atom diameter is %f\n", pPixel->GetMeanCoreAtomDiameter() );
113
114		// Precompute alphaOne
115		double alphaOne = 1 - eigenvalues[0];
116
117		// Watch out for /0 problems
118		if(alphaOne==0)
119		{
120		alphaOne = 0.001;
121		}
122		// Now compute the rest of the lengths
123		for(unsigned int i = 1; i < NDimensions; i++)
124		{
125	LEN,SEM	axisLengthArray[i] = ( (1 - eigenvalues[i]) / alphaOne) * axisLengthArray[0] ;
126		}
127
128		// Build the ellipsoid voting region
129	LEN	typedef EllipsoidInteriorExteriorSpatialFunction<NDimensions, PositionType> VoteFunctionType;
130		typename VoteFunctionType::Pointer ellipsoid = VoteFunctionType::New();
131
132		// Create an iterator to traverse the ellipsoid region
133		typedef FloodFilledSpatialFunctionConditionalIterator
134	IND	******<Self, VoteFunctionType> ItType;
135
136	LEN	// The seed position for the ellipsoid is the current pixel's index in data space
137		// since this is always at the center of the voting ellipsoid
138		typename Self::IndexType seedPos = bloxIt.GetIndex();
139
140		// Figure out the center of the ellipsoid, which is the center
141		// of the voting pixel
142		typename VoteFunctionType::InputType centerPosition;
143
144		ContinuousIndex<double, dim> contIndex;
145
146		for(unsigned int i = 0; i < dim; i ++ )
147		{
148		contIndex[i] = (double)seedPos[i] + 0.5;
149		}
150
151		// Get the physical location of this center index
152		this->TransformContinuousIndexToPhysicalPoint(contIndex, centerPosition);
153
154		ellipsoid->SetCenter(centerPosition);
155		ellipsoid->SetOrientations(eigenvectors);
156		ellipsoid->SetAxes(axisLengthArray);
157
158		// Instantiate the iterator
159		ItType sfi = ItType(this, ellipsoid, seedPos);
160
161		// Get the position of the voting blox
162		typedef Point<double, NDimensions> TPosition;
163		TPosition voterPosition;
164		typename Self::IndexType voterIndex = bloxIt.GetIndex();
165		this->TransformIndexToPhysicalPoint(voterIndex, voterPosition);
166
167		int voteeCount = 0;
168
169		sfi.SetCenterInclusionStrategy();
170
171		// Iterate through the ellipsoid and cast votes
172		for( sfi.GoToBegin(); !( sfi.IsAtEnd() ); ++sfi)
173		{
174		TPosition voteePosition;
175		typename Self::IndexType voteeIndex = sfi.GetIndex();
176		//std::cerr << "voteeIndex "<< voteeIndex << std::endl ;
177
178		this->TransformIndexToPhysicalPoint(voteeIndex, voteePosition);
179
180		// vector from voting blox to current votee
181		typename TPosition::VectorType dbar = voterPosition - voteePosition;
182		voteeCount ++;
183
184		// The voting process and variables are explained in
185		// IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 18, NO. 10, OCTOBER 1999
186		// page 1029
187
188		// The votee does not get voted for if it's empty
189		if( sfi.Get().GetSize() == 0 )
190		{
191		continue;
192		}
193		// form the ellipsoidal distance de
194		double de = 0;
195		double sf_de_sqr = 0;
196
197		for (unsigned int i = 0; i < NDimensions; i++)
198		{
199	LEN	de += pow((dot_product(eigenvectors.get_column(i), dbar.GetVnlVector() ) /
200		axisLengthArray[i] ), 2);
201		}
202
203		de = sqrt(de);
204
205		//printf("De = %f\n", de);
206
207		double weight_factor = exp(-1.0*de*de);
208
209		// vote strength
210		double voteStrength = pPixel->size()*weight_factor;
211		//printf("Vote strength = %f\n", voteStrength);
212		//printf("weight_factor = %f\n", weight_factor);
213
214		// Get eigenanalysis results
215		sf_eigenvalues = sfi.Get().GetEigenvalues();
216		sf_eigenvectors = sfi.Get().GetEigenvectors();
217
218		for (unsigned int i = 0; i < NDimensions; i++)
219		{
220	LEN	sf_de_sqr += pow((dot_product(sf_eigenvectors.get_column(i), dbar.GetVnlVector() ) /
221		axisLengthArray[i] ), 2);
222		}
223
224		//printf("sf_de = %f\n", sqrt(sf_de_sqr));
225
226		//CALCULATE WEIGHT FACTOR FOR INDEX OF SPATIAL FUNCTION ITERATION
227		double sf_weight_factor = exp(-1.0*sf_de_sqr);
228
229	LEN	//HERE WE CALL CalcWeightedCoreAtomLocation using de to keep track of the weighted location of
230		//each voted medial node based on constituent core atom locations
231		// Reminder: sfi.Get() is the pixel being voted for
232		// and pPixel is doing the voting
233		sfi.Get().CalcWeightedCoreAtomLocation(sf_weight_factor, pPixel);
234
235		// cast the vote
236	LEN	sfi.Get().CollectVote(pPixel->GetRawCMatrixPointer(), voteStrength, pPixel->size() );
237
238		} // end cast votes from this pixel
239
240		//printf("Blox voted for %i other pixels\n", voteeCount);
241		} // end cast votes from all pixels
242
243		// The final task is to normalize all of the voted blox
244		// and recompute the eigenanalysis on the new matrix
245		for(bloxIt.GoToBegin(); !bloxIt.IsAtEnd(); ++bloxIt)
246		{
247		(&bloxIt.Value())->NormalizeVotedCMatrix();
248		(&bloxIt.Value())->DoVotedEigenanalysis();
249		}
250
251		m_MedialNodeCount = voterCount;
252		//cerr << "MedialNodeCount = " << m_MedialNodeCount << endl;
253		}
254
255		template <unsigned int dim>
256		void
257		BloxCoreAtomImage<dim>
258		::PrintSelf(std::ostream& os, Indent indent) const
259		{
260		Superclass::PrintSelf(os,indent);
261
262		// Iterator to access all pixels in the image
263		ImageRegionConstIterator<Self> bloxIt =
264	IND	****ImageRegionConstIterator<Self>(this, this->GetLargestPossibleRegion() );
265
266		// Pointer for accessing pixels
267		BloxCoreAtomPixel<NDimensions> pPixel;
268
269		// Results of eigenanalysis from each pixel
270		typename BloxCoreAtomPixel<NDimensions>::EigenvalueType eigenvalues;
271		typename BloxCoreAtomPixel<NDimensions>::EigenvalueType veigenvalues;
272
273	LEN	os << indent << "Index\t# Core Atoms\tEigen Values\t\t\tMean CA Length\tVoted Eigen Values\n"
274	LEN	<< "-----\t------------\t------------\t\t\t--------------\t------------------\n" << std::endl;
275
276		int counter = 0;
277
278		for(bloxIt.GoToBegin(); !bloxIt.IsAtEnd(); ++bloxIt)
279		{
280		// Get a pointer to the pixel
281		pPixel = bloxIt.Value();
282
283		eigenvalues = pPixel.GetEigenvalues();
284		veigenvalues = pPixel.GetVotedEigenvalues();
285
286		if(!pPixel.empty())
287		{
288		os << indent << bloxIt.GetIndex() << "\t";
289		os << indent << pPixel.GetSize() << "\t";
290	LEN	os << indent << eigenvalues[0] << " " << eigenvalues[1] << " " << eigenvalues[2] << "\t";
291		os << indent << pPixel.GetMeanCoreAtomDiameter() << "\t\t";
292	LEN	os << indent << veigenvalues[0] << " " << veigenvalues[1] << " " << veigenvalues[2] << "\t" << std::endl;
293	LEN	os << std::endl << indent << "Node Pointer List: " << (*m_NodePointerList)[counter]->GetVotedLocation() << std::endl;
294		counter++;
295		}
296		}
297		os << "Number of Medial Nodes: " << m_MedialNodeCount << std::endl;
298		}
299
300		} // end namespace itk
301
302		#endif
303

---

Generated by KWStyle 1.0b on Tuesday January,17 at 02:14:00PM

© Kitware Inc.