Commit b07ee870 authored by aylward's avatar aylward
Browse files

ENH: Direct access to tutorials from within Slicer3



git-svn-id: http://svn.slicer.org/Slicer4/trunk@4342 3bd1e089-480b-0410-8dfb-8563597acbee
parent c9d60a81
...@@ -370,6 +370,15 @@ TARGET_LINK_LIBRARIES (${CLP}Lib ITKIO ITKBasicFilters) ...@@ -370,6 +370,15 @@ TARGET_LINK_LIBRARIES (${CLP}Lib ITKIO ITKBasicFilters)
SUBDIRS ( Realign ) SUBDIRS ( Realign )
SUBDIRS( DiffusionApplications ) SUBDIRS( DiffusionApplications )
SUBDIRS( ExtractSkeleton )
IF(USE_MIDAS)
SUBDIRS( MIDASApplications )
ENDIF(USE_MIDAS)
#IF(USE_BatchMake)
#SUBDIRS( BatchMakeApplications )
#ENDIF(USE_BatchMake)
IF (Slicer3_SOURCE_DIR) IF (Slicer3_SOURCE_DIR)
# install each target in the production area (where it would appear in an # install each target in the production area (where it would appear in an
......
PROJECT (ExtractSkeleton)
# Disable MSVC 8 warnings
IF(WIN32)
OPTION(DISABLE_MSVC8_DEPRECATED_WARNINGS
"Disable Visual Studio 8 deprecated warnings" ON)
MARK_AS_ADVANCED(FORCE DISABLE_MSVC8_DEPRECATED_WARNINGS)
IF(DISABLE_MSVC8_DEPRECATED_WARNINGS)
ADD_DEFINITIONS(-D_CRT_SECURE_NO_DEPRECATE)
ENDIF(DISABLE_MSVC8_DEPRECATED_WARNINGS)
ENDIF(WIN32)
SET ( ExtractSkeleton_SOURCE
ExtractSkeleton.cxx
SkelGraph.h
SkelGraph.cxx
tilg_iso_3D.h
tilg_iso_3D.cxx
coordTypes.h
misc.h
misc.cxx )
GENERATECLP(ExtractSkeleton_SOURCE ExtractSkeleton.xml)
ADD_EXECUTABLE(ExtractSkeleton ${ExtractSkeleton_SOURCE})
TARGET_LINK_LIBRARIES(ExtractSkeleton ITKIO)
ADD_LIBRARY(ExtractSkeletonModule SHARED ${ExtractSkeleton_SOURCE})
SET_TARGET_PROPERTIES (ExtractSkeletonModule PROPERTIES
COMPILE_FLAGS "-Dmain=ModuleEntryPoint")
TARGET_LINK_LIBRARIES(ExtractSkeletonModule ITKIO)
IF (Slicer3_SOURCE_DIR)
SET(TARGET_MODULES ExtractSkeletonModule ExtractSkeleton)
FOREACH(target ${TARGET_MODULES})
INSTALL(TARGETS ${target}
RUNTIME DESTINATION ${SLICER_INSTALL_LIBRARIES_DIR}/Plugins
LIBRARY DESTINATION ${SLICER_INSTALL_LIBRARIES_DIR}/Plugins)
ENDFOREACH(target)
ENDIF (Slicer3_SOURCE_DIR)
/*=========================================================================
Program: Extract Skeleton
Module: $RCSfile: main.cxx,v $
Language: C++
Date: $Date: 2007-07-10 11:35:36 -0400 (Tue, 10 Jul 2007) $
Version: $Revision: 48 $
Copyright (c) 2007 Insight Software Consortium. All rights reserved.
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
#ifdef __BORLANDC__
#define ITK_LEAN_AND_MEAN
#endif
#include <cmath>
#include <iostream>
#include "itkImage.h"
#include "itkImageFileWriter.h"
#include "itkImageFileReader.h"
#include "itkImageRegionIterator.h"
#include "itkImageRegionConstIterator.h"
#include "SkelGraph.h"
#include "coordTypes.h"
#include "tilg_iso_3D.h"
#include "ExtractSkeletonCLP.h"
/** Main command */
int main(int argc, char **argv)
{
std::cout << "Begun." << std::endl;
PARSE_ARGS;
std::cout << "Processed args." << std::endl;
try
{
typedef unsigned char InputPixelType;
typedef itk::Image<InputPixelType,3> InputImageType;
typedef itk::ImageFileReader<InputImageType> ReaderType;
typedef unsigned char OutputPixelType;
typedef itk::Image<OutputPixelType, 3> OutputImageType;
typedef itk::ImageFileWriter<OutputImageType> WriterType;
ReaderType::Pointer Reader = ReaderType::New();
Reader->SetFileName(InputImageFileName.c_str());
Reader->Update();
InputImageType::Pointer inputImage = Reader->GetOutput();
std::cout << "Read image." << std::endl;
OutputImageType::Pointer outputImage = OutputImageType::New();
outputImage->SetRegions(Reader->GetOutput()->GetLargestPossibleRegion());
outputImage->SetSpacing(Reader->GetOutput()->GetSpacing());
outputImage->SetOrigin(Reader->GetOutput()->GetOrigin());
outputImage->Allocate();
itk::Size<3> itkSize = inputImage->GetLargestPossibleRegion().GetSize();
int dim [3];
dim[0] = itkSize[0];
dim[1] = itkSize[1];
dim[2] = itkSize[2];
InputPixelType * inputImageBuffer = inputImage->GetBufferPointer();
OutputPixelType * outputImageBuffer = outputImage->GetBufferPointer();
memset(outputImageBuffer, 0, dim[0]*dim[1]*dim[2]*sizeof(OutputPixelType));
std::cout << "Initialized output image." << std::endl;
int extract2DSheet = 0;
if(SkeletonType == "2D")
{
extract2DSheet = 1;
}
tilg_iso_3D(dim[0], dim[1], dim[2],
inputImageBuffer, outputImageBuffer, extract2DSheet);
std::cout << "Extracted skeleton." << std::endl;
std::list<point> axisPoints;
SkelGraph * graph = new SkelGraph();
graph->Extract_skel_graph(outputImageBuffer, dim);
graph->Extract_max_axis_in_graph();
graph->Sample_along_axis(NumberOfPoints, &axisPoints);
std::ofstream writeOutputFile;
writeOutputFile.open(OutputPointsFileName.c_str());
if(!DontPruneBranches)
{
memset(outputImageBuffer, 0,
dim[0]*dim[1]*dim[2]*sizeof(OutputPixelType));
}
int i = 0;
std::list<point>::iterator iter = axisPoints.begin();
while(iter != axisPoints.end())
{
OutputImageType::IndexType pt;
pt[0] = iter->x;
pt[1] = iter->y;
pt[2] = iter->z;
if(!DontPruneBranches)
{
outputImage->SetPixel(pt, 255);
}
writeOutputFile << i << " " << iter->x
<< " " << iter->y
<< " " << iter->z << std::endl;
iter++;
i++;
}
std::cout << "Wrote points file." << std::endl;
WriterType::Pointer writer = WriterType::New();
writer->SetFileName(OutputImageFileName.c_str());
writer->SetInput(outputImage);
writer->Update();
std::cout << "Wrote output image." << std::endl;
}
catch( itk::ExceptionObject &excep )
{
std::cerr << argv[0] << " : Exception caught!" << std::endl;
std::cerr << excep << std::endl;
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
<?xml version="1.0" encoding="utf-8"?>
<executable>
<category>Filtering</category>
<title>Extract Skeleton</title>
<description>
Extract the skeleton of a binary object. The skeleton can be limited to being a 1D curve or allowed to be a full 2D manifold. The branches of the skeleton can be pruned so that only the maximal center skeleton is returned.
</description>
<version>0.1.0.$Revision: 2104 $(alpha)</version>
<documentation-url></documentation-url>
<license></license>
<contributor>Pierre Seroul, Martin Styner, Guido Gerig, and Stephen Aylward
</contributor>
<acknowledgements>
The original implementation of this method was provided by ETH Zurich, Image Analysis Laboratory of Profs Olaf Kuebler, Gabor Szekely and Guido Gerig. Martin Styner at UNC, Chapel Hill made enhancements. Wrapping for Slicer was provided by Pierre Seroul and Stephen Aylward at Kitware, Inc.
</acknowledgements>
<parameters>
<label>IO</label>
<description>Input/output parameters</description>
<image>
<name>InputImageFileName</name>
<label>Input Image</label>
<channel>input</channel>
<index>1</index>
<description>Input image</description>
</image>
<image>
<name>OutputImageFileName</name>
<label>Output Image</label>
<channel>output</channel>
<index>2</index>
<description>Skeleton of the input image</description>
</image>
</parameters>
<parameters>
<label>Skeleton</label>
<description>Skeleton parameters</description>
<string-enumeration>
<name>SkeletonType</name>
<longflag>type</longflag>
<label>Skeleton type</label>
<description>Type of skeleton to create</description>
<default>1D</default>
<element>1D</element>
<element>2D</element>
</string-enumeration>
<boolean>
<name>DontPruneBranches</name>
<longflag>dontPrune</longflag>
<label>Do not prune branches</label>
<description>Return the full skeleton, not just the maximal skeleton</description>
<default>false</default>
</boolean>
<integer>
<name>NumberOfPoints</name>
<longflag>numPoints</longflag>
<label>Number Of Points</label>
<description>Number of points used to represent the skeleton</description>
<default>100</default>
</integer>
<string>
<name>OutputPointsFileName</name>
<longflag>pointsFile</longflag>
<label>Output points file (.txt)</label>
<description>Name of the file to store the coordinates of the central (1D) skeleton points
</description>
<default>skeleton.txt</default>
</string>
</parameters>
</executable>
This diff is collapsed.
#ifndef _SKEL_GRAPH_H_
#define _SKEL_GRAPH_H_
#include <list>
#include "coordTypes.h"
typedef struct skel_branch_struct
{
int branchID; // == position in graph
double length;
double acc_length; // for temporary use when searching maximal path
list<int> * acc_path;
double max_length;
list<int> * max_path; // maximal path
point * end_1_point;
point * end_2_point;
list<int> * end_1_neighbors; // id's == one can use advance for random access
list<int> * end_2_neighbors;
} skel_branch;
class SkelGraph {
private:
// Extracted Graph
list<skel_branch> * graph;
// To_Do list, only of temporary use for extract graph
list<skel_branch> * to_do;
// endpoint list, only of temporary use
list<point> * endpoints;
// Image to extract from
unsigned char *image;
int dim[3];
// Label image, only of temporary use
int *label_image;
skel_branch * max_node; // for storage of start of maximal path
double max_length;
// private routines
void Add_new_elem_to_todo(skel_branch * &newElem);
// adds a new element with default values to To_do list
void find_endpoints();
// find all endpoints in image
int endpoint_Test(int x, int y, int z);
// tests whether (x,y,z) is an endpoint
void get_valid_neighbors(point *point1, std::list<point> * &neighbors);
// returns a list of valid neighbors at act_point
void ResetGraph();
public:
SkelGraph();
SkelGraph(SkelGraph * graph); // not fully implemented
~SkelGraph();
void PrintGraph();
void Extract_skel_graph(unsigned char *orig_image, int orig_dim[3]);
// Extract skeletal graph
void Extract_max_axis_in_graph();
// extract maximal path between 2 points in the graph
void Sample_along_axis(int n_dim, list<point> * axis_points);
// sample along the medial axis and perpendicular to it
};
#endif
/*****************************************************************************************************
//
// Title: coordTypes.h
//
*******************************************************************************************************/
#ifndef COORD_TYPES_H
#define COORD_TYPES_H
#include "math.h"
#include <vector>
#include "misc.h"
using namespace std;
typedef struct point_struct
{
int x;
int y;
int z;
} point;
class Coord3i
{
int v[3];
public:
inline int& operator[](int i) { return v[i]; }
inline void conv(double * i) {i[0] = v[0]; i[1] = v[1]; i[2] = v[2];};
};
class Coord3f
{
float v[3];
public:
inline float& operator[](int i) { return v[i]; }
inline void conv(float * i) {i[0] = v[0]; i[1] = v[1]; i[2] = v[2];};
inline void conv(double * i) {i[0] = v[0]; i[1] = v[1]; i[2] = v[2];};
};
class Coord3d
{
double v[3];
public:
inline double& operator[](int i) { return v[i]; };
inline void conv(int * i) {i[0] = (int) v[0]; i[1] = (int) v[1]; i[2] = (int) v[2];};
inline void conv(float * i) {i[0] = v[0]; i[1] = v[1]; i[2] = v[2];};
inline void conv(double * i) {i[0] = v[0]; i[1] = v[1]; i[2] = v[2];};
};
inline void normcrossprod(double *v1, double *v2, double *norm)
// calculate normalized crossproduct
{
double absval;
norm[0] = v1[1] * v2[2] - v1[2] * v2[1];
norm[1] = v1[2] * v2[0] - v1[0] * v2[2];
norm[2] = v1[0] * v2[1] - v1[1] * v2[0];
absval = sqrt(norm[0]*norm[0] + norm[1]*norm[1] + norm[2]*norm[2]);
norm[0] /= absval;
norm[1] /= absval;
norm[2] /= absval;
}
inline double vectorangle(double *v1, double *v2)
// calculate angle between two vectors (0..M_PI), you might want to adjust to 0..M_PI/2
// range after call via 'if (angle > M_PI/2) angle = M_PI-angle;'
{
double prod = 0, length1 = 0, length2 = 0;
for (int k = 0; k < 3; k++) {
prod += v1[k]*v2[k];
length1 += v1[k]*v1[k];
length2 += v2[k]*v2[k];
}
return acos(prod/sqrt(length1*length2));
}
inline double vectorangle(Coord3d v1, Coord3d v2)
// calculate angle between two vectors (0..M_PI), you might want to adjust to 0..M_PI/2
// range after call via 'if (angle > M_PI/2) angle = M_PI-angle;'
{
double prod = 0, length1 = 0, length2 = 0;
for (int k = 0; k < 3; k++) {
prod += v1[k]*v2[k];
length1 += v1[k]*v1[k];
length2 += v2[k]*v2[k];
}
return acos(prod/sqrt(length1*length2));
}
inline double vec_length(Coord3d x) {
return sqrt(sqr(x[0]) + sqr(x[1]) + sqr(x[2]));
}
inline double vec_length(double *x) {
return sqrt(sqr(x[0]) + sqr(x[1]) + sqr(x[2]));
}
inline double vec_length(double *x, double *y) {
return sqrt(sqr(x[0]-y[0]) + sqr(x[1]-y[1]) + sqr(x[2]-y[2]));
}
inline int transWorldToImage(Coord3d loc_world, int *loc_img,
double *origin, int *dims, double voxelsize)
// transform and check index, returns 0 on success and 1 if corrected location
{
int adjust = 0;
for (int i = 0; i < 3; i++) {
loc_img [i] = (int) ((loc_world[i] - origin[i]) / voxelsize);
if (loc_img[i] < 0 ) { adjust = 1; loc_img[i] = 0; }
if (loc_img[i] >= dims[i] ) { loc_img[i] = dims[i] - 1;adjust = 1; }
}
return adjust;
}
inline int transWorldToImage(double *loc_world, int *loc_img,
double *origin, int *dims, double voxelsize)
// transform and check index, returns 0 on success and 1 if corrected location
{
int adjust = 0;
for (int i = 0; i < 3; i++) {
loc_img [i] = (int) ((loc_world[i] - origin[i]) / voxelsize);
if (loc_img[i] < 0 ) { adjust = 1; loc_img[i] = 0; }
if (loc_img[i] >= dims[i] ) { loc_img[i] = dims[i] - 1;adjust = 1; }
}
return adjust;
}
#endif
/*
* misc.cc -- miscellaneous utitlity functions
*
* author: msturm
* created: 27 Mar 1997
* changes: mastyner
*/
#include "misc.h"
void* ipAllocateData(const int size, const size_t elemsize)
{
void *data = NULL;
if (!(data = malloc(size * elemsize))) {
fprintf(stderr,"Error: ipAllocateData [%s, line %d]: memory allocation failed:",
__FILE__, __LINE__);
perror("");
exit(errno);
}
memset(data, 0, size * elemsize);
return(data);
}
size_t ipGetDataSize(const ipDataType type)
{
size_t retval;
switch (type) {
case IP_BYTE:
retval = sizeof(char);
break;
case IP_SHORT:
retval = sizeof(short);
break;
case IP_INT:
retval = sizeof(int);
break;
case IP_FLOAT:
retval = sizeof(float);
break;
case IP_DOUBLE:
retval = sizeof(double);
break;
default:
#ifdef DEBUG_VSKEL
fprintf(stderr,"Warning: ipGetDataSize [%s, line %d]: unsuported data type: %d\n",
__FILE__, __LINE__, type);
#endif
retval = 0;
break;
}
return retval;
}
/*
* misc.h
*
* author: msturm
* created: 27 Mar 1997
* changes: mastyner
*/
#ifndef __IP_MISC_H__
#define __IP_MISC_H__
#include <stdlib.h>
#include <stdio.h>
#include <errno.h>
#include <sys/types.h>
#include <string.h>
typedef enum {
IP_BYTE = 0, /* AVS_TYPE_BYTE = 0 */
IP_INT, /* AVS_TYPE_INTEGER = 1 */
IP_FLOAT, /* AVS_TYPE_REAL = 2 */
IP_DOUBLE, /* AVS_TYPE_DOUBLE = 3 */
IP_SHORT /* AVS_TYPE_SHORT = 4 */
} ipDataType;