root/Testing/TutorialTest.cxx

Revision 1, 18.3 kB (checked in by bob, 6 years ago)

Importing ITKGray

Line 
1 /*=========================================================================
2
3   Program:   Insight Segmentation & Registration Toolkit
4   Module:    $RCSfile: TutorialTest.cxx,v $
5   Language:  C++
6   Date:      $Date: 2006/12/02 04:22:20 $
7   Version:   $Revision: 1.1 $
8   Copyright (c) 2003 Insight Consortium. All rights reserved.
9   See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
10
11      This software is distributed WITHOUT ANY WARRANTY; without even
12      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
13      PURPOSE.  See the above copyright notices for more information.
14 =========================================================================*/
15 // Borland compiler is very lazy so we need to instantiate the template
16 //  by hand
17 #if defined(__BORLANDC__)
18 #include <SNAPBorlandDummyTypes.h>
19 #endif
20
21 #include "ImageIORoutines.h"
22 #include "IRISApplication.h"
23 #include "IRISImageData.h"
24 #include "LevelSetMeshPipeline.h"
25 #include "SNAPImageData.h"
26
27 #include "vtkPolyData.h"
28 #include "vtkPolyDataWriter.h"
29
30 #include <vector>
31
32 /**
33  * This function is designed to serve as a tutorial through the classes
34  * that form the logical framework of SNAP. In this test, the student will
35  * be able to load images, perform preprocessing, initialize and run segmentation,
36  * output slices from the segmentation and relabel the segmentatinon
37  *
38  * The method should be invoked with two parameters. First is the full path to
39  * the file MRIcrop-orig.gipl which contains the test image. Second is the full
40  * path to the label file MRIcrop-seg.label. Both of these files can be downloaded
41  * from the SNAP website
42  */
43 int main(int argc, char *argv[])
44 {
45   if(argc < 3)
46     {
47     cerr << "Must specify two arguments." << endl;
48     return 1;
49     }
50   /* =======================================================================
51    * Section 1. Starting the Application and Loading an Image
52    * ====================================================================== */
53
54   /* We begin by creating an IRISApplication object. This is the most high
55    * level object in the SNAP class hierarchy and it provides access to the
56    * remainder of the SNAP logic classes */
57   IRISApplication application;
58
59   /* Next we will load the test image. In this example it is expected that the
60    * image MRIcrop-orig.gipl, which can be downloaded from the SNAP website
61    * will be used. The user will pass the full path to this file as the first
62    * parameter to the program */
63   const char *sImageFileName = argv[1];
64  
65   try {
66     /* The responsibility for loading the image is with the user. First, we
67      * create an itk::Image object into which the image will be loaded from
68      * disk */
69     IRISApplication::GreyImageType::Pointer imgGrey;
70
71     /* SNAP provides a helper method LoadImageFromFile for loading images */
72     LoadImageFromFile(sImageFileName,imgGrey);
73
74     /* Once the image has been loaded, it can be passed on to the IRISApplication.
75      * However, SNAP needs to know the orientation of the image, which is
76      * expressed as the position of the image's origin {X,Y,Z}={0,0,0} in
77      * patient coordinates. This information is passed on in the form of a
78      * three letter RAI code. The value "ASR" means that the origin lies in the
79      * Anterior-Superior-Right corner of the patient coordinate system, with
80      * X axis corresponding to the Anterior-Posterior axis, Y axis corresponding
81      * to the Inferior-Superior axis and Z to the Right-Left axis. For the image
82      * we are loading, the correct code is "RAI" */
83     application.UpdateIRISGreyImage(imgGrey,"RAI");
84
85     /* In addition to loading the image, we are going to load the segmentation
86      * labels associated with this image. Segmentation labels are stored in the
87      * file MRIcrop-seg.label and passed on as an argument to this test program */
88     const char *sLabelFileName = argv[2];
89     application.GetColorLabelTable()->LoadFromFile(sLabelFileName);
90   }
91   catch(itk::ExceptionObject &exception)
92   {
93     // This code is called if LoadImageFromFile fails.
94     cout << "Error Loading Image: " << exception;
95     return -1;
96   }
97
98   /* =======================================================================
99    * Section 2. Accessing Image Data Using IRISImageData Class
100    * ====================================================================== */
101                                      
102   /* Now the image has been loaded into IRISApplication. IRISApplication stores
103    * its data using two high level classes: IRISImageData and SNAPImageData.
104    * IRISImageData holds the original greyscale image that we have loaded as
105    * well as the labelling of that image that is constructed throughout the
106    * segmentation. SNAPImageData is used for automatic segmentation and contains
107    * a subregion of the original greyscale image. At this point SNAPImageData has
108    * not yet been allocated. SNAPImageData is a child class of IRISImageData.
109    *
110    * The method IRISApplication::GetCurrentImageData() returns the IRISImageData
111    * object associated with the current state of SNAP. In manual state (current
112    * state right now) it returns a pointer to an IRISImageData object. In
113    * automatic segmentation state, it returns a pointer to a SNAPImageData
114    * object.
115    *
116    * The IRISImageData and SNAPImageData objects can also be accessed explicitly
117    * using the GetIRISImageData() and GetSNAPImageData() pointers.
118    *
119    * We will begin by examining the IRISImageData class. We can use it to access
120    * the underlying image and to examine the segmentation. */
121   IRISImageData *irisData = application.GetCurrentImageData();
122
123   /* IRISImageData contains two images: the grey level image and a label or
124    * segmentation image. Both images are of the same dimensions. The grey image
125    * has voxels of type GreyType and the segmentation image of type LabelType */
126   cout << "Image Dimensions: " <<
127     to_int(irisData->GetVolumeExtents()) << endl;       // Ignore the to_int!
128 //  cout << "Voxel Size      : " << irisData->GetVoxelScaleFactor() << endl;
129
130   /* IRISImageData also handles segmentation labels. We can inquire about each label.
131    * Labels are indexed from 1 to 255 (0 is the clear label). Only some labels are
132    * defined as 'valid', i.e., available to the user to edit */
133   for(unsigned int iLabel=0; iLabel < MAX_COLOR_LABELS; iLabel++)
134     {
135     /* Labels are represented by the class ColorLabel. Each label has a color,
136      * a transparency level, a visible flag, a 'visible in 3d mesh flag', and a
137      * textual description */
138     const ColorLabel &label = application.GetColorLabelTable()->GetColorLabel(iLabel);                                                       
139     cout << "Label Number  : " << iLabel << endl;
140     cout << "  Decription  : " << label.GetLabel() << endl;
141     cout << "  Color       : {"
142       << (int) label.GetRGB(0) << ","
143       << (int) label.GetRGB(1) << ","
144       << (int) label.GetRGB(2) << "}" << endl;
145     cout << "  Opacity     : " << label.GetAlpha() << endl;
146     }
147
148   /* The current drawing label and the current draw-over label are properties of
149    * the SNAP application and are accessed using the GlobalState class.
150    * GlobalState stores a great number of application-wide settings */
151   application.GetGlobalState()->SetDrawingColorLabel(7);  // caudate label
152   application.GetGlobalState()->SetCoverageMode(PAINT_OVER_ALL);
153
154   /* The next class we will learn about is the ImageWrapper class, which is a
155    * wrapper around the standard itk::Image class. The ImageWrapper provides
156    * several mini-pipelines involving the itk::Image stored within. First
157    * let's access the grey images' wrapper in IRISImageData */
158   GreyImageWrapper *wrapGrey = irisData->GetGrey();
159
160   /* The class GreyImageWrapper is a subclass of the generic templated class
161    * ImageWrapper<TPixel> with template argument TPixel=GreyType. We can
162    * use ImageWrapper to get to the itk::Image itself */
163   GreyImageWrapper::ImagePointer imgGrey = wrapGrey->GetImage();
164
165   /* Here are some of the methods available through the ImageWrapper class */
166   cout << "Grey Min        : " << wrapGrey->GetImageMin() << endl;
167   cout << "Grey Max        : " << wrapGrey->GetImageMax() << endl;
168
169   /* The IRISImageData and ImageWrapper classes keep track of the current
170    * position of the SNAP 3D cursor. When an image is loaded, the cursor is
171    * positioned in the center of the image, as we can see by querying the
172    * ImageWrapper class */
173   Vector3ui vCursor = wrapGrey->GetSliceIndex();
174   cout << "Cursor Position : " << to_int(vCursor) << endl; // Ignore the to_int!
175   cout << "Value at Cursor : "
176     << wrapGrey->GetVoxel(vCursor) << endl;
177
178   /* The cursor position obtained above is given in image coordinates, where
179    * Z is the slice number, Y is the row number and  X is the column number.
180    * In order to specify slice coordinates in patient (anatomical)
181    * coordinates, we need to know the correct transformation, which we had \
182    * specified using an 'RAI' code when loading the grey image. The
183    * transformation can be accessed through the IRISApplication object */
184   ImageCoordinateGeometry geometry = irisData->GetImageGeometry();
185   Vector3ui vCursorPatient =
186     geometry.GetImageToAnatomyTransform().TransformVoxelIndex(vCursor);
187   cout << "Anatomical Posn.: " << to_int(vCursorPatient) << endl;
188
189   /* The cursor position is the same for the grey image wrapper and the label
190    * image wrapper. Hence, to set the cursor position, we must call the
191    * appropriate method in IRISImageData. We can also call the method
192    * IRISApplication::SetCursorPosition for the same effect. */
193   vCursor -= Vector3ui(1,0,0);
194   irisData->SetCrosshairs(vCursor);
195
196   /* ImageWrapper provides a way to extract orthogonal 2D slices form the
197    * image. A slice is extracted by setting the crosshair position and calling
198    * the GetDisplaySlice method. This returns a 2D image of type unsigned char
199    * that can be displayed on the screen.
200    *
201    * The mapping between the intensities in the grey image (which are shorts)
202    * and intensities in the slice (which are bytes) can be changed by calling
203    * the method SetIntensityMapFunction() */
204   GreyImageWrapper::DisplaySlicePointer imgSlice = wrapGrey->GetDisplaySlice(0);
205
206   /* So far, we have looked at the grey image wrapper. The segmentation image
207    * wrapper is very similar, except that it produces slices that are RGB
208    * images.
209    *
210    * When we loaded the grey image, the segmentation image in IRISImageData was
211    * automatically initialized to the same size and filled with voxels of
212    * label 0, which is the 'Clear' label in SNAP. To update the segmentation,
213    * we can directly modify the image wrapped by the LabelImageWrapper or
214    * we can use the GetVoxelForUpdate() method: */
215   irisData->GetSegmentation()->GetVoxelForUpdate(vCursor) = 2;
216   LabelImageWrapper::DisplaySlicePointer imgRGBSlice =
217     irisData->GetSegmentation()->GetDisplaySlice(0);
218
219   /* For validation purposes, we are going to write to disk the grey and
220    * segmentation slices that we've extracted in this Section */
221   SaveImageToFile("greyslice01.png",imgSlice.GetPointer());
222   SaveImageToFile("rgbslice01.png",imgRGBSlice.GetPointer());
223  
224   /* =======================================================================
225    * Section 3. Automatic Segmentation Mode
226    * ====================================================================== */
227                                      
228   /* To begin automatic segmentation, we need to switch on the automatic
229    * segmentation mode. In order to do that, we must specify the region of
230    * interest that we want to pass on to this mode and whether or not we
231    * want to perform resampling on the region of interest. In this case we
232    * will use the entire image as the region of interest and do no resampling */
233   SNAPSegmentationROISettings roiSettings;
234   roiSettings.SetROI(irisData->GetImageRegion());
235   roiSettings.SetResampleFlag(false);
236  
237   /* We can now tell the IRISApplication to initialize the SNAPImageData object.
238    * This method takes the ROI settings as the first parameter and an
239    * itk::Command object as the optional second parameter. The Command is used
240    * as a callback during the resampling, and is provided primarify for GUIs
241    * that want to display a progress bar. */
242   application.InitializeSNAPImageData(roiSettings);
243
244   /* Once the SNAP data is initalized, we can instruct the IRISApplication to
245    * use it as the current image data, in other words, to switch to the auto
246    * segmentation mode */
247   application.SetCurrentImageDataToSNAP();
248
249   /* The SNAPImageData class inherits the methods and attributes of the
250    * IRISImageData class and provides additional functionality related to
251    * level set image segmentation. In addition to the Grey and Segmentation
252    * images (and image wrappers), the SNAPImageData class contains a Speed
253    * image wrapper, a SnakeInitialization image wrapper and a Snake image
254    * wrapper. All three of these image wrappers contain itk::Image objects
255    * with pixel type float. In addition, the SNAPImageData class provides
256    * method for running the automatic segmetnation pipeline */
257   SNAPImageData *snapData = application.GetSNAPImageData();
258
259   /* Automatic segmentation begins with preprocessing. We must choose the type
260    * of preprocessing to apply and the parameters of the preprocessing. In this
261    * tutorial we will use edge-based snakes and edge-detection preprocessing */
262   EdgePreprocessingSettings edgeSettings =
263     EdgePreprocessingSettings::MakeDefaultSettings();
264   edgeSettings.SetGaussianBlurScale(0.6);
265   edgeSettings.SetRemappingSteepness(0.02);
266   edgeSettings.SetRemappingExponent(2.0);
267  
268   /* The preprocessing is performed when we call the DoEdgePreprocessing method.
269    * the first argument is an EdgePreprocessingSettings object and the second
270    * optional argument is a Command that can be used to implement a progress bar*/
271   snapData->DoEdgePreprocessing(edgeSettings);
272
273   /* The result of the preprocessing is stored in the Speed image wrapper in
274    * the SNAPImageData class. We can write out a slice of the preprocessed image */
275   SpeedImageWrapper *wrapSpeed = snapData->GetSpeed();
276   SaveImageToFile("speedslice00.mha",wrapSpeed->GetDisplaySlice(0));
277
278   /* The next step after preprocessing is to initialize the segmentation with
279    * bubbles. This is performed by passing an array of bubbles to the method
280    * SNAPImageData::InitializeSegmentationPipeline. The bubbles given below
281    * are located in the caudate nuclei of our image */
282   std::vector<Bubble> bubbles(2);
283   bubbles[0].center = Vector3i(50,29,26); bubbles[0].radius = 3;
284   bubbles[1].center = Vector3i(50,57,34); bubbles[1].radius = 2;
285
286   /* In addition to the bubble array, the method InitializeSegmentationPipeline
287    * requires a set of level set segmentation parameters, which are stored in
288    * the SnakeParameters class */
289   SnakeParameters parameters = SnakeParameters::GetDefaultEdgeParameters();
290
291   /* In order for the snake to converge, we enable the advection force */
292   parameters.SetAdvectionWeight(5.0);
293
294   /* We can now initialize the segmentation pipeline. In addition to the
295    * snake parameters and the bubbles, we need to specify the color label to
296    * be used for the snake-based segmentation */
297   snapData->InitializeSegmentation(
298     parameters,bubbles,application.GetGlobalState()->GetDrawingColorLabel());
299
300   /* Now the pipeline is initialized and ready to run. Run 250 iterations */
301   snapData->RunSegmentation(250);
302
303   /* We can also rewind the segmentation */
304   snapData->RestartSegmentation();
305
306   /* And we can change the parameters on the fly */
307   snapData->RunSegmentation(100);
308   parameters.SetAdvectionWeight(3.0);
309   snapData->SetSegmentationParameters(parameters);
310   snapData->RunSegmentation(200);
311
312   /* Now the segmentation is done. The Snake image wrapper contains a floating
313    * point image whose positive voxels correspond to the pixels outside of the
314    * segmentation boundary and whose negative valued pixels are inside. In
315    * to get an actual snake surface, we can use the LevelSetMeshPipeline object,
316    * which uses VTK to trace the zero-level-set of the Snake image */ 
317   LevelSetMeshPipeline meshPipeline;   
318   meshPipeline.SetImage(snapData->GetSnake()->GetImage());
319   // meshPipeline.SetImage(snapData->GetLevelSetImage());
320
321   /* Pass the globally stored mesh options to the mesh pipeline */
322   meshPipeline.SetMeshOptions(application.GetGlobalState()->GetMeshOptions());
323  
324   /* Create a vtkPolyData to store the resulting mesh, and run the pipeline */
325   vtkPolyData *meshResult = vtkPolyData::New(); 
326   meshPipeline.ComputeMesh(meshResult);
327
328   /* Export the mesh to VTK format */
329   vtkPolyDataWriter *polyWriter = vtkPolyDataWriter::New();
330   polyWriter->SetInput(meshResult);
331   polyWriter->SetFileTypeToASCII();
332   polyWriter->SetFileName("levelset.poly");
333   polyWriter->Write();
334   polyWriter->Delete();
335
336   /* Great. We have now executed automatic segmentation and exported the leve set
337    * as a mesh. What's left is to exit the automatic segmentation mode and to
338    * merge the segmentation results with the multi-label segmentation image
339    * stored in IRISImageData. The following method will do that for us. */
340   application.UpdateIRISWithSnapImageData();
341   application.SetCurrentImageDataToIRIS();
342
343   /* Since we are done with the segmentation data, we should release resources
344    * associated with it */
345   application.ReleaseSNAPImageData();
346
347   /* =======================================================================
348    * Section 4. Additional Functionality
349    * ====================================================================== */
350                                      
351   /* One of the features of SNAP it to be able to divide segmentations into two
352    * labels using a cut-plane. This can be done programatically, using the
353    * following code. In our example, this will divide the left and the right
354    * caudates, giving them different labels */
355
356   /* We will paint with label 9 over label 8 */
357   application.GetGlobalState()->SetDrawingColorLabel(9); 
358   application.GetGlobalState()->SetCoverageMode(PAINT_OVER_COLORS);
359   application.GetGlobalState()->SetOverWriteColorLabel(8); 
360
361   /* This actually performs the cut, given an normal vector and intercept that
362    * define the cut plane */
363   application.RelabelSegmentationWithCutPlane(Vector3d(1,0,0), 64);
364
365   /* We can now save the segmentation result as a 3d image */
366   SaveImageToFile("result.gipl",irisData->GetSegmentation()->GetImage());
367
368   /* Thank you for reading this tutorial! */
369   return 0;
370 }
371
372 // A global variable required for linkage
373 std::ostream &verbose = std::cout;
Note: See TracBrowser for help on using the browser.