Commit d47497b8 authored by Ken Martin's avatar Ken Martin
Browse files

some fixes and cleaning

parent 05e7f55c
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
Date: $Date$ Date: $Date$
Version: $Revision$ Version: $Revision$
Copyright (c) 1993-1995 Ken Martin, Will Schroeder, Bill Lorensen. Copyright (c) 1993-1997 Ken Martin, Will Schroeder, Bill Lorensen.
This software is copyrighted by Ken Martin, Will Schroeder and Bill Lorensen. This software is copyrighted by Ken Martin, Will Schroeder and Bill Lorensen.
The following terms apply to all files associated with the software unless The following terms apply to all files associated with the software unless
...@@ -37,107 +37,51 @@ MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. ...@@ -37,107 +37,51 @@ MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
=========================================================================*/ =========================================================================*/
#include <math.h> #include <stdlib.h>
#include "vtkLinkEdgels.h" #include "vtkLinkEdgels.h"
#include "vtkMath.h"
//----------------------------------------------------------------------------
// Description: // Description:
// Construct instance of vtkImageLinkEdgels with GradientThreshold set to // Construct instance of vtkLinkEdgels with GradientThreshold set to
// 0.1, PhiThreshold set to 90 degrees and LinkThreshold set to 90 degrees. // 0.1, PhiThreshold set to 90 degrees and LinkThreshold set to 90 degrees.
vtkLinkEdgels::vtkLinkEdgels() vtkLinkEdgels::vtkLinkEdgels()
{ {
this->Input = NULL;
this->GradientThreshold = 0.1; this->GradientThreshold = 0.1;
this->PhiThreshold = 90; this->PhiThreshold = 90;
this->LinkThreshold = 90; this->LinkThreshold = 90;
} }
//----------------------------------------------------------------------------
// Description:
// This filter executes if it or a previous filter has been modified or
// if its data has been released and it is forced to update.
void vtkLinkEdgels::Update()
{
int execute;
// make sure input is available
if ( !this->Input )
{
vtkErrorMacro(<< "No input...can't execute!");
return;
}
execute = this->Input->GetPipelineMTime() > this->ExecuteTime
|| this->GetMTime() > this->ExecuteTime
|| this->Output->GetDataReleased();
if (execute)
{
vtkDebugMacro(<< "ConditionalUpdate: Condition satisfied, executeTime = "
<< this->ExecuteTime
<< ", modifiedTime = " << this->GetMTime()
<< ", input MTime = " << this->Input->GetPipelineMTime()
<< ", released = " << this->Output->GetDataReleased());
if ( this->StartMethod ) (*this->StartMethod)(this->StartMethodArg);
this->Output->Initialize(); //clear output
this->Execute();
this->ExecuteTime.Modified();
this->SetDataReleased(0);
if ( this->EndMethod ) (*this->EndMethod)(this->EndMethodArg);
}
}
void vtkLinkEdgels::Execute() void vtkLinkEdgels::Execute()
{ {
vtkImageRegion *region = new vtkImageRegion; vtkPointData *pd;
int regionExtent[8];
vtkFloatPoints *newPts=0; vtkFloatPoints *newPts=0;
vtkCellArray *newLines=0; vtkCellArray *newLines=0;
vtkFloatScalars *inScalars;
vtkFloatScalars *outScalars; vtkFloatScalars *outScalars;
vtkStructuredPoints *input = this->GetInput();
int numPts;
vtkFloatVectors *outVectors; vtkFloatVectors *outVectors;
vtkPolyData *output = this->GetOutput(); vtkPolyData *output = this->GetOutput();
int sliceNum, sliceMin, sliceMax; int *dimensions;
float *CurrMap, *inDataPtr;
// error checking vtkVectors *inVectors;
if ( this->Input == NULL ) int ptId;
{
vtkErrorMacro(<<"Execute:Please specify an input!");
return;
}
// Fill in image information.
this->Input->UpdateImageInformation(region);
region->SetAxes(VTK_IMAGE_X_AXIS, VTK_IMAGE_Y_AXIS,
VTK_IMAGE_COMPONENT_AXIS, VTK_IMAGE_Z_AXIS);
// get the input region vtkDebugMacro(<< "Extracting structured points geometry");
region->GetImageExtent(4, regionExtent);
region->SetExtent(4, regionExtent);
this->Input->UpdateRegion(region);
if ( ! region->AreScalarsAllocated()) pd = input->GetPointData();
dimensions = input->GetDimensions();
inScalars = (vtkFloatScalars *)pd->GetScalars();
inVectors = pd->GetVectors();
if ((numPts=this->Input->GetNumberOfPoints()) < 2 || inScalars == NULL)
{ {
vtkErrorMacro(<< "Execute: Could not get region."); vtkErrorMacro(<<"No data to transform!");
return; return;
} }
// chekc data type for float // set up the input
if (region->GetScalarType() != VTK_FLOAT) inDataPtr = inScalars->GetPtr(0);
{
vtkImageRegion *temp = region;
vtkWarningMacro(<<"Converting non float image data to float");
region = new vtkImageRegion;
region->SetScalarType(VTK_FLOAT);
region->SetExtent(VTK_IMAGE_DIMENSIONS, temp->GetExtent());
region->CopyRegionData(temp);
temp->Delete();
}
// Finally do edge following to extract the edge data from the Thin image // Finally do edge following to extract the edge data from the Thin image
newPts = new vtkFloatPoints; newPts = new vtkFloatPoints;
newLines = new vtkCellArray; newLines = new vtkCellArray;
...@@ -146,22 +90,15 @@ void vtkLinkEdgels::Execute() ...@@ -146,22 +90,15 @@ void vtkLinkEdgels::Execute()
vtkDebugMacro("doing edge linking\n"); vtkDebugMacro("doing edge linking\n");
// //
// for each slice link edgels // Traverse all points, for each point find Gradient in the Image map.
// //
sliceMin = regionExtent[6]; for (ptId=0; ptId < dimensions[2]; ptId++)
sliceMax = regionExtent[7];
for (sliceNum = sliceMin; sliceNum <= sliceMax; sliceNum++)
{ {
// Set extent to be just one slice. CurrMap = inDataPtr + dimensions[0]*dimensions[1]*ptId;
regionExtent[6] = regionExtent[7] = sliceNum;
region->SetExtent(4, regionExtent); this->LinkEdgels(dimensions[0],dimensions[1],CurrMap, inVectors,
this->LinkEdgels(region, newLines,newPts,outScalars,outVectors, sliceNum); newLines,newPts,outScalars,outVectors,ptId);
} }
// restore original extent.
regionExtent[6] = sliceMin;
regionExtent[7] = sliceMax;
region->SetExtent(4, regionExtent);
output->SetPoints(newPts); output->SetPoints(newPts);
output->SetLines(newLines); output->SetLines(newLines);
...@@ -175,28 +112,25 @@ void vtkLinkEdgels::Execute() ...@@ -175,28 +112,25 @@ void vtkLinkEdgels::Execute()
newLines->Delete(); newLines->Delete();
outScalars->Delete(); outScalars->Delete();
outVectors->Delete(); outVectors->Delete();
region->Delete();
} }
// Description: // Description:
// This method links the edges for one image. // This method links the edges for one image.
void vtkLinkEdgels::LinkEdgels(vtkImageRegion *region, void vtkLinkEdgels::LinkEdgels(int xdim, int ydim, float *image,
vtkCellArray *newLines, vtkVectors *inVectors,
vtkFloatPoints *newPts, vtkCellArray *newLines,
vtkFloatScalars *outScalars, vtkFloatPoints *newPts,
vtkFloatVectors *outVectors, vtkFloatScalars *outScalars,
int z) vtkFloatVectors *outVectors,
int z)
{ {
int *extent = region->GetExtent();
int xdim = extent[1] - extent[0] + 1;
int ydim = extent[3] - extent[2] + 1;
int **forward; int **forward;
int **backward; int **backward;
int x,y,ypos,zpos; int x,y,ypos,zpos;
int currX, currY, i; int currX, currY, i;
int newX, newY; int newX, newY;
int startX, startY; int startX, startY;
float vec[3], vec1[2], vec2[2]; float vec[3], vec1[3], vec2[3];
float linkThresh, phiThresh; float linkThresh, phiThresh;
// these direction vectors are rotated 90 degrees // these direction vectors are rotated 90 degrees
// to convert gradient direction into edgel direction // to convert gradient direction into edgel direction
...@@ -208,11 +142,9 @@ void vtkLinkEdgels::LinkEdgels(vtkImageRegion *region, ...@@ -208,11 +142,9 @@ void vtkLinkEdgels::LinkEdgels(vtkImageRegion *region,
static int xoffset[8] = {1,1,0,-1,-1,-1,0,1}; static int xoffset[8] = {1,1,0,-1,-1,-1,0,1};
static int yoffset[8] = {0,1,1,1,0,-1,-1,-1}; static int yoffset[8] = {0,1,1,1,0,-1,-1,-1};
int length, start; int length, start;
int bestDirection = 0; int bestDirection;
float error, bestError; float error, bestError;
float *imgPtrX, *imgPtrY, *imgPtrX2;
int imgIncX, imgIncY, imgIncVec;
forward = new int *[ydim]; forward = new int *[ydim];
backward = new int *[ydim]; backward = new int *[ydim];
for (i = 0; i < ydim; i++) for (i = 0; i < ydim; i++)
...@@ -227,19 +159,15 @@ void vtkLinkEdgels::LinkEdgels(vtkImageRegion *region, ...@@ -227,19 +159,15 @@ void vtkLinkEdgels::LinkEdgels(vtkImageRegion *region,
linkThresh = cos(this->LinkThreshold*3.1415926/180.0); linkThresh = cos(this->LinkThreshold*3.1415926/180.0);
phiThresh = cos(this->PhiThreshold*3.1415926/180.0); phiThresh = cos(this->PhiThreshold*3.1415926/180.0);
imgPtrY = (float *)region->GetScalarPointer();
region->GetIncrements(imgIncX,imgIncY,imgIncVec);
// first find all forward & backwards links // first find all forward & backwards links
for (y = 0; y < ydim; y++, imgPtrY += imgIncY) for (y = 0; y < ydim; y++)
{ {
ypos = y*xdim; ypos = y*xdim;
imgPtrX = imgPtrY; for (x = 0; x < xdim; x++)
for (x = 0; x < xdim; x++, imgPtrX += imgIncX)
{ {
// find forward and backward neighbor for this pixel // find forward and backward neighbor for this pixel
// if its value is less than threshold then ignore it // if its value is less than threshold then ignore it
if ((*imgPtrX) < this->GradientThreshold) if (image[x+ypos] < this->GradientThreshold)
{ {
forward[y][x] = -1; forward[y][x] = -1;
backward[y][x] = -1; backward[y][x] = -1;
...@@ -247,8 +175,8 @@ void vtkLinkEdgels::LinkEdgels(vtkImageRegion *region, ...@@ -247,8 +175,8 @@ void vtkLinkEdgels::LinkEdgels(vtkImageRegion *region,
else else
{ {
// try all neighbors as forward, first try four connected // try all neighbors as forward, first try four connected
vec1[0] = *(imgPtrX + imgIncVec); inVectors->GetVector(x+ypos+zpos,vec1);
vec1[1] = *(imgPtrX + 2*imgIncVec); vtkMath::Normalize(vec1);
// first eliminate based on phi1 - alpha // first eliminate based on phi1 - alpha
bestError = 0; bestError = 0;
for (i = 0; i < 8; i += 2) for (i = 0; i < 8; i += 2)
...@@ -262,13 +190,13 @@ void vtkLinkEdgels::LinkEdgels(vtkImageRegion *region, ...@@ -262,13 +190,13 @@ void vtkLinkEdgels::LinkEdgels(vtkImageRegion *region,
if ((x + xoffset[i] >= 0)&&(x + xoffset[i] < xdim)&& if ((x + xoffset[i] >= 0)&&(x + xoffset[i] < xdim)&&
(y + yoffset[i] >= 0)&&(y + yoffset[i] < ydim)&& (y + yoffset[i] >= 0)&&(y + yoffset[i] < ydim)&&
(!backward[y+yoffset[i]][x+xoffset[i]])&& (!backward[y+yoffset[i]][x+xoffset[i]])&&
(*(imgPtrX + xoffset[i]*imgIncX + yoffset[i]*imgIncY) >= (image[x + xoffset[i] + (y+yoffset[i])*xdim] >=
this->GradientThreshold)) this->GradientThreshold))
{ {
// satisfied the first test, now check second // satisfied the first test, now check second
imgPtrX2 = imgPtrX + xoffset[i]*imgIncX + yoffset[i]*imgIncY; inVectors->GetVector(x + xoffset[i] +
vec2[0] = *(imgPtrX2 + imgIncVec); (y + yoffset[i])*xdim + zpos,vec2);
vec2[1] = *(imgPtrX2 + 2*imgIncVec); vtkMath::Normalize(vec2);
if ((vec1[0]*vec2[0] + vec1[1]*vec2[1]) >= phiThresh) if ((vec1[0]*vec2[0] + vec1[1]*vec2[1]) >= phiThresh)
{ {
// pased phi - phi test does the forward neighbor // pased phi - phi test does the forward neighbor
...@@ -310,13 +238,13 @@ void vtkLinkEdgels::LinkEdgels(vtkImageRegion *region, ...@@ -310,13 +238,13 @@ void vtkLinkEdgels::LinkEdgels(vtkImageRegion *region,
if ((x + xoffset[i] >= 0)&&(x + xoffset[i] < xdim)&& if ((x + xoffset[i] >= 0)&&(x + xoffset[i] < xdim)&&
(y + yoffset[i] >= 0)&&(y + yoffset[i] < ydim)&& (y + yoffset[i] >= 0)&&(y + yoffset[i] < ydim)&&
(!backward[y+yoffset[i]][x+xoffset[i]])&& (!backward[y+yoffset[i]][x+xoffset[i]])&&
(*(imgPtrX + xoffset[i]*imgIncX + yoffset[i]*imgIncY) >= (image[x + xoffset[i] + (y+yoffset[i])*xdim] >=
this->GradientThreshold)) this->GradientThreshold))
{ {
// satisfied the first test, now check second // satisfied the first test, now check second
imgPtrX2 = imgPtrX + xoffset[i]*imgIncX + yoffset[i]*imgIncY; inVectors->GetVector(x + xoffset[i] +
vec2[0] = *(imgPtrX2 + imgIncVec); (y + yoffset[i])*xdim + zpos,vec2);
vec2[1] = *(imgPtrX2 + 2*imgIncVec); vtkMath::Normalize(vec2);
if ((vec1[0]*vec2[0] + vec1[1]*vec2[1]) >= phiThresh) if ((vec1[0]*vec2[0] + vec1[1]*vec2[1]) >= phiThresh)
{ {
// pased phi - phi test does the forward neighbor // pased phi - phi test does the forward neighbor
...@@ -351,8 +279,6 @@ void vtkLinkEdgels::LinkEdgels(vtkImageRegion *region, ...@@ -351,8 +279,6 @@ void vtkLinkEdgels::LinkEdgels(vtkImageRegion *region,
// now construct the chains // now construct the chains
imgPtrX = (float *)region->GetScalarPointer();
vec[2] = z; vec[2] = z;
for (y = 0; y < ydim; y++) for (y = 0; y < ydim; y++)
{ {
...@@ -385,12 +311,10 @@ void vtkLinkEdgels::LinkEdgels(vtkImageRegion *region, ...@@ -385,12 +311,10 @@ void vtkLinkEdgels::LinkEdgels(vtkImageRegion *region,
{ {
currX = newX; currX = newX;
currY = newY; currY = newY;
outScalars->InsertNextScalar(*(imgPtrX + currX*imgIncX + outScalars->InsertNextScalar(image[currX + currY*xdim]);
currY*imgIncY)); inVectors->GetVector(currX+currY*xdim+zpos,vec2);
vtkMath::Normalize(vec2);
vec[0] = *(imgPtrX + currX*imgIncX + currY*imgIncY + imgIncVec); outVectors->InsertNextVector(vec2);
vec[1] = *(imgPtrX + currX*imgIncX + currY*imgIncY + 2*imgIncVec);
outVectors->InsertNextVector(vec);
vec[0] = currX; vec[0] = currX;
vec[1] = currY; vec[1] = currY;
newPts->InsertNextPoint(vec); newPts->InsertNextPoint(vec);
...@@ -431,10 +355,9 @@ void vtkLinkEdgels::LinkEdgels(vtkImageRegion *region, ...@@ -431,10 +355,9 @@ void vtkLinkEdgels::LinkEdgels(vtkImageRegion *region,
void vtkLinkEdgels::PrintSelf(ostream& os, vtkIndent indent) void vtkLinkEdgels::PrintSelf(ostream& os, vtkIndent indent)
{ {
vtkPolySource::PrintSelf(os,indent); vtkStructuredPointsToPolyDataFilter::PrintSelf(os,indent);
os << indent << "GradientThreshold:" << this->GradientThreshold << "\n"; os << indent << "GradientThreshold:" << this->GradientThreshold << "\n";
os << indent << "LinkThreshold:" << this->LinkThreshold << "\n"; os << indent << "LinkThreshold:" << this->LinkThreshold << "\n";
os << indent << "PhiThreshold:" << this->PhiThreshold << "\n"; os << indent << "PhiThreshold:" << this->PhiThreshold << "\n";
} }
...@@ -6,7 +6,6 @@ ...@@ -6,7 +6,6 @@
Date: $Date$ Date: $Date$
Version: $Revision$ Version: $Revision$
Copyright (c) 1993-1995 Ken Martin, Will Schroeder, Bill Lorensen. Copyright (c) 1993-1995 Ken Martin, Will Schroeder, Bill Lorensen.
This software is copyrighted by Ken Martin, Will Schroeder and Bill Lorensen. This software is copyrighted by Ken Martin, Will Schroeder and Bill Lorensen.
...@@ -38,27 +37,48 @@ MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. ...@@ -38,27 +37,48 @@ MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
=========================================================================*/ =========================================================================*/
// .NAME vtkLinkeEdgles - takes a gradient image and links edgels // .NAME vtkLinkEdgels - links edgels together to form digital curves.
// .SECTION Description
// vtkLinkEdgels links edgels into digital curves which are then stored
// as polylines. The algorithm works one pixel at a time only looking at
// its immediate neighbors. There is a GradientThreshold that can be set
// that eliminates any pixels with a smaller gradient value. This can
// be used as the lower threshold of a two value edgel thresholding.
//
// For the remaining edgels, links are first tried for the four
// connected neighbors. A succesful neighbor will satisfy three
// tests. First both edgels must be above the gradient
// threshold. Second, the difference between the orientation between
// the two edgels (Alpha) and each edgels orientation (Phi) must be
// less than LinkThreshold. Third, the difference between the two
// edgels Phi values must be less than PhiThreshold.
// The most successful link is selected. The meaure is simply the
// sum of the three angle differences (actually stored as the sum of
// the cosines). If none of the four connect neighbors succeds, then
// the eight connect neighbors are examined using the same method.
//
// This filter requires gradient information so you will need to use
// a vtkImageGradient at some point prior to this filter. Typically
// a vtkNonMaximumSuppression filter is also used. vtkThresholdEdgels
// can be used to complete the two value edgel thresholding as used
// in a Canny edge detector. The vtkSubpixelPositionEdgels filter
// can also be used after this filter to adjust the edgel locations.
// .SECTION see also
// vtkImage vtkImageGradient vtkNonMaximumSuppression
#ifndef __vtkLinkEdgels_h #ifndef __vtkLinkEdgels_h
#define __vtkLinkEdgels_h #define __vtkLinkEdgels_h
#include "vtkPolySource.h" #include "vtkStructuredPointsToPolyDataFilter.h"
#include "vtkImageSource.h"
#include "vtkImageRegion.h"
class vtkLinkEdgels : public vtkPolySource class vtkLinkEdgels : public vtkStructuredPointsToPolyDataFilter
{ {
public: public:
vtkLinkEdgels(); vtkLinkEdgels();
char *GetClassName() {return "vtkImageToStructurePoints";}; char *GetClassName() {return "vtkLinkEdgels";};
void PrintSelf(ostream& os, vtkIndent indent); void PrintSelf(ostream& os, vtkIndent indent);
// Description:
// Set/Get the input object from the image pipline.
vtkSetObjectMacro(Input,vtkImageSource);
vtkGetObjectMacro(Input,vtkImageSource);
// Description: // Description:
// Set/Get the threshold for Phi vs. Alpha link thresholding. // Set/Get the threshold for Phi vs. Alpha link thresholding.
vtkSetMacro(LinkThreshold,float); vtkSetMacro(LinkThreshold,float);
...@@ -74,21 +94,15 @@ public: ...@@ -74,21 +94,15 @@ public:
vtkSetMacro(GradientThreshold,float); vtkSetMacro(GradientThreshold,float);
vtkGetMacro(GradientThreshold,float); vtkGetMacro(GradientThreshold,float);
void Update();
protected: protected:
vtkImageSource *Input; void Execute();
void LinkEdgels(int xdim, int ydim,float *image, vtkVectors *inVectors,
vtkCellArray *newLines, vtkFloatPoints *newPts,
vtkFloatScalars *outScalars, vtkFloatVectors *outVectors,
int z);
float GradientThreshold; float GradientThreshold;
float PhiThreshold; float PhiThreshold;
float LinkThreshold; float LinkThreshold;
void LinkEdgels(vtkImageRegion *region,
vtkCellArray *newLines, vtkFloatPoints *newPts,
vtkFloatScalars *outScalars,
vtkFloatVectors *outVectors, int z);
void Execute();
}; };
#endif #endif
<
...@@ -7,7 +7,7 @@ ...@@ -7,7 +7,7 @@
Version: $Revision$ Version: $Revision$
Copyright (c) 1993-1995 Ken Martin, Will Schroeder, Bill Lorensen. Copyright (c) 1993-1996 Ken Martin, Will Schroeder, Bill Lorensen.
This software is copyrighted by Ken Martin, Will Schroeder and Bill Lorensen. This software is copyrighted by Ken Martin, Will Schroeder and Bill Lorensen.
The following terms apply to all files associated with the software unless The following terms apply to all files associated with the software unless
...@@ -48,97 +48,53 @@ MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. ...@@ -48,97 +48,53 @@ MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
// 0.1, PhiThreshold set to 90 degrees and LinkThreshold set to 90 degrees. // 0.1, PhiThreshold set to 90 degrees and LinkThreshold set to 90 degrees.
vtkLinkSurfels::vtkLinkSurfels() vtkLinkSurfels::vtkLinkSurfels()
{ {
this->Input = NULL;
this->GradientThreshold = 0.1; this->GradientThreshold = 0.1;
this->PhiThreshold = 90; this->PhiThreshold = 90;
this->LinkThreshold = 90; this->LinkThreshold = 90;
} }
//----------------------------------------------------------------------------
// Description:
// This filter executes if it or a previous filter has been modified or
// if its data has been released and it is forced to update.
void vtkLinkSurfels::Update()
{
int execute;
// make sure input is available
if ( !this->Input )
{
vtkErrorMacro(<< "No input...can't execute!");
return;
}
execute = this->Input->GetPipelineMTime() > this->ExecuteTime
|| this->GetMTime() > this->ExecuteTime
|| this->Output->GetDataReleased();
if (execute)
{
vtkDebugMacro(<< "Update: Condition satisfied, executeTime = "
<< this->ExecuteTime
<< ", modifiedTime = " << this->GetMTime()
<< ", input MTime = " << this->Input->GetPipelineMTime()
<< ", released = " << this->Output->GetDataReleased());
if ( this->StartMethod ) (*this->StartMethod)(this->StartMethodArg);
this->Output->Initialize(); //clear output
this->Execute();
this->ExecuteTime.Modified();
this->SetDataReleased(0);
if ( this->EndMethod ) (*this->EndMethod)(this->EndMethodArg);
}
}
void vtkLinkSurfels::Execute() void vtkLinkSurfels::Execute()
{ {
vtkImageRegion *region = new vtkImageRegion; vtkPointData *pd;
int regionBounds[8];
vtkFloatPoints *newPts=0; vtkFloatPoints *newPts=0;