/*========================================================================= Program: Visualization Toolkit Module: $RCSfile: vtkMyImplicitModeller.cxx,v $ Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen All rights reserved. See Copyright.txt or http://www.kitware.com/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 notice for more information. =========================================================================*/ #include "vtkMyImplicitModeller.h" #include "vtkCell.h" #include "vtkCellLocator.h" #include "vtkClipPolyData.h" #include "vtkCommand.h" #include "vtkFloatArray.h" #include "vtkGenericCell.h" #include "vtkImageData.h" #include "vtkInformation.h" #include "vtkInformationVector.h" #include "vtkMath.h" #include "vtkMultiThreader.h" #include "vtkMutexLock.h" #include "vtkObjectFactory.h" #include "vtkPlane.h" #include "vtkPointData.h" #include "vtkPolyData.h" #include "vtkRectilinearGrid.h" #include "vtkStreamingDemandDrivenPipeline.h" #include "vtkStructuredGrid.h" #include "vtkUnstructuredGrid.h" #include "vtkImageIterator.h" #include "vtkImageProgressIterator.h" #include vtkCxxRevisionMacro(vtkMyImplicitModeller, "$Revision: 1.98 $"); vtkStandardNewMacro(vtkMyImplicitModeller); struct vtkMyImplicitModellerAppendInfo { vtkMyImplicitModeller *Modeller; vtkDataSet **Input; double MaximumDistance; }; //---------------------------------------------------------------------------- // Construct with sample dimensions=(50,50,50), and so that model bounds are // automatically computed from the input. Capping is turned on with CapValue // equal to a large positive number. vtkMyImplicitModeller::vtkMyImplicitModeller() { this->MaximumDistance = 0.1; this->ModelBounds[0] = 0.0; this->ModelBounds[1] = 0.0; this->ModelBounds[2] = 0.0; this->ModelBounds[3] = 0.0; this->ModelBounds[4] = 0.0; this->ModelBounds[5] = 0.0; this->BoundsComputed = 0; this->SampleDimensions[0] = 50; this->SampleDimensions[1] = 50; this->SampleDimensions[2] = 50; this->Capping = 1; this->OutputScalarType = VTK_FLOAT; this->CapValue = this->GetScalarTypeMax( this->OutputScalarType ); this->ScaleToMaximumDistance = 0; // only used for non-float output type this->DataAppended = 0; this->AdjustBounds = 1; this->AdjustDistance = 0.0125; this->ProcessMode = VTK_CELL_MODE; this->LocatorMaxLevel = 5; this->Threader = vtkMultiThreader::New(); this->NumberOfThreads = this->Threader->GetNumberOfThreads(); } vtkMyImplicitModeller::~vtkMyImplicitModeller() { if (this->Threader) { this->Threader->Delete(); } } void vtkMyImplicitModeller::SetOutputScalarType(int type) { double scalarMax; vtkDebugMacro(<< this->GetClassName() << " (" << this << "): setting OutputScalarType to " << type); scalarMax = this->GetScalarTypeMax(type); if (scalarMax) // legal type { int modified = 0; if (this->CapValue != scalarMax) { this->CapValue = scalarMax; modified = 1; } if (this->OutputScalarType != type) { this->OutputScalarType = type; modified = 1; } if (modified) { this->Modified(); } } } void vtkMyImplicitModeller::SetCapValue(double value) { vtkDebugMacro(<< this->GetClassName() << " (" << this << "): setting CapValue to " << value); // clamp to between 0 and max for scalar type double max = this->GetScalarTypeMax(this->OutputScalarType); if (this->CapValue != (value < 0 ? 0 : (value > max ? max : value))) { this->CapValue = (value < 0 ? 0 : (value > max ? max : value)); this->Modified(); } } double vtkMyImplicitModeller::GetScalarTypeMax(int type) { switch (type) { case VTK_UNSIGNED_CHAR: return (double)VTK_UNSIGNED_CHAR_MAX; case VTK_CHAR: return (double)VTK_CHAR_MAX; case VTK_UNSIGNED_SHORT: return (double)VTK_UNSIGNED_SHORT_MAX; case VTK_SHORT: return (double)VTK_SHORT_MAX; case VTK_UNSIGNED_INT: return (double)VTK_UNSIGNED_INT_MAX; case VTK_INT: return (double)VTK_INT_MAX; case VTK_UNSIGNED_LONG: return (double)VTK_UNSIGNED_LONG_MAX; case VTK_LONG: return (double)VTK_LONG_MAX; case VTK_FLOAT: return (double)VTK_FLOAT_MAX; case VTK_DOUBLE: return (double)VTK_DOUBLE_MAX; default: return 0; } } //---------------------------------------------------------------------------- // Initialize the filter for appending data. You must invoke the // StartAppend() method before doing successive Appends(). It's also a // good idea to manually specify the model bounds; otherwise the input // bounds for the data will be used. void vtkMyImplicitModeller::StartAppend() { this->StartAppend(0); } void vtkMyImplicitModeller::StartAppend(int internal) { vtkIdType numPts; vtkIdType i; double maxDistance; if (!internal) { // we must call update information because we can't be sure that // it has been called. this->UpdateInformation(); } this->GetOutput()->SetUpdateExtent(this->GetOutput()->GetWholeExtent()); vtkDebugMacro(<< "Initializing data"); this->AllocateOutputData(this->GetOutput()); this->UpdateProgress(0.0); this->DataAppended = 1; numPts = this->SampleDimensions[0] * this->SampleDimensions[1] * this->SampleDimensions[2]; // initialize output to CapValue at each location maxDistance = this->CapValue; vtkDataArray *newScalars = this->GetOutput()->GetPointData()->GetScalars(); for (i=0; iSetComponent(i, 0, maxDistance); } } template void SetOutputDistance(double distance, OT *outputValue, double capValue, double scaleFactor) { // for now, just doing "normal" cast... could consider doing round? if (scaleFactor) // need to scale the distance { *outputValue = static_cast(distance * scaleFactor); } else { if (capValue && distance > capValue) // clamping iff non-float type { distance = capValue; } *outputValue = static_cast(distance); } } // Convert distance as stored in output (could be scaled and/or non-double // type) to double distance with correct scaling void ConvertToDoubleDistance(double inDistance, double &distance, double &distance2, double scaleFactor) { if (scaleFactor) { distance = inDistance * scaleFactor; } else { distance = inDistance; } distance2 = distance * distance; } //---------------------------------------------------------------------------- // Templated append for VTK_VOXEL_MODE process mode and any type of output data template void vtkMyImplicitModellerAppendExecute(vtkMyImplicitModeller *self, vtkDataSet *input, vtkImageData *outData, int outExt[6], double maxDistance, vtkCellLocator *locator, int id, OT *) { int i, j, k; int subId; vtkIdType cellId; double pcoords[3]; double *spacing, *origin; double maxDistance2 = maxDistance * maxDistance; double x[3], prevDistance, prevDistance2, distance2, betterDistance; double closestPoint[3], mDist; // allocate weights for the EvaluatePosition double *weights = new double[input->GetMaxCellSize()]; // Traverse each voxel; using CellLocator to find the closest point vtkGenericCell *cell = vtkGenericCell::New(); spacing = outData->GetSpacing(); origin = outData->GetOrigin(); vtkImageProgressIterator outIt(outData, outExt, self, id); // so we know how to scale if desired double scaleFactor = 0; // 0 used to indicate not scaling double toDoubleScaleFactor = 0; // 0 used to indicate not scaling double capValue = 0; // 0 used to indicate not clamping (float or double) if (self->GetOutputScalarType() != VTK_FLOAT && self->GetOutputScalarType() != VTK_DOUBLE) { capValue = self->GetCapValue(); if (self->GetScaleToMaximumDistance()) { scaleFactor = capValue / maxDistance; toDoubleScaleFactor = maxDistance / capValue; } } int testIndex = 0; for (k = outExt[4]; k <= outExt[5]; k++) { x[2] = spacing[2] * k + origin[2]; for (j = outExt[2]; j <= outExt[3]; j++) { cellId = -1; x[1] = spacing[1] * j + origin[1]; OT* outSI = outIt.BeginSpan(); for (i = outExt[0]; i <= outExt[1]; i++, testIndex++) { x[0] = spacing[0] * i + origin[0]; ConvertToDoubleDistance(*outSI, prevDistance, prevDistance2, toDoubleScaleFactor); betterDistance = -1; if (cellId != -1) { cell->EvaluatePosition(x, closestPoint, subId, pcoords, distance2, weights); if (distance2 <= maxDistance2 && distance2 < prevDistance2) { mDist = sqrt(distance2); betterDistance = mDist; } else if (prevDistance2 < maxDistance2) { mDist = prevDistance; } else { mDist = maxDistance; } } else if (prevDistance2 < maxDistance2) { mDist = prevDistance; } else { mDist = maxDistance; } if (locator->FindClosestPointWithinRadius(x, mDist, closestPoint, cell, cellId, subId, distance2) ) { if(distance2 <= prevDistance2) { betterDistance = sqrt(distance2); } } else { cellId = -1; } if (betterDistance != -1) { SetOutputDistance(betterDistance, outSI, capValue, scaleFactor); } outSI++; } outIt.NextSpan(); } } cell->Delete(); delete [] weights; } //---------------------------------------------------------------------------- // This is the multithreaded piece of the append when doing per voxel // processing - it is called once for each thread, with each thread // taking a different slab of the output to work on. The acutal work is done // in vtkMyImplicitModellerAppendExecute; here we just setup for the per voxel // processing. static VTK_THREAD_RETURN_TYPE vtkMyImplicitModeller_ThreadedAppend( void *arg ) { int threadCount; int threadId; vtkMyImplicitModellerAppendInfo *userData; vtkImageData *output; double maxDistance; int i; double *bounds, adjBounds[6]; double *spacing; double *origin; int slabSize, slabMin, slabMax; int outExt[6]; threadId = ((vtkMultiThreader::ThreadInfo *)(arg))->ThreadID; threadCount = ((vtkMultiThreader::ThreadInfo *)(arg))->NumberOfThreads; userData = (vtkMyImplicitModellerAppendInfo *) (((vtkMultiThreader::ThreadInfo *)(arg))->UserData); if (userData->Input[threadId] == NULL) { return VTK_THREAD_RETURN_VALUE; } maxDistance = userData->MaximumDistance; output = userData->Modeller->GetOutput(); spacing = output->GetSpacing(); origin = output->GetOrigin(); int *sampleDimensions = userData->Modeller->GetSampleDimensions(); if (!output->GetPointData()->GetScalars()) { vtkGenericWarningMacro("Sanity check failed."); return VTK_THREAD_RETURN_VALUE; } // break up into slabs based on threadId and threadCount slabSize = sampleDimensions[2] / threadCount; if (slabSize == 0) // in case threadCount > sampleDimensions[2] { slabSize = 1; } slabMin = threadId * slabSize; if (slabMin >= sampleDimensions[2]) { return VTK_THREAD_RETURN_VALUE; } slabMax = slabMin + slabSize - 1; if (threadId == threadCount - 1) { slabMax = sampleDimensions[2] - 1; } bounds = userData->Input[threadId]->GetBounds(); for (i=0; i<3; i++) { adjBounds[2*i] = bounds[2*i] - maxDistance; adjBounds[2*i+1] = bounds[2*i+1] + maxDistance; } // compute dimensional bounds in data set for (i = 0; i < 3; i++) { outExt[i*2] = (int) ((double)(adjBounds[2*i] - origin[i]) / spacing[i]); outExt[i*2+1] = (int) ((double)(adjBounds[2*i+1] - origin[i]) / spacing[i]); if (outExt[i*2] < 0) { outExt[i*2] = 0; } if (outExt[i*2+1] >= sampleDimensions[i]) { outExt[i*2+1] = sampleDimensions[i] - 1; } } // input not close enough to effect this slab if (outExt[4] > slabMax || outExt[5] < slabMin) { return VTK_THREAD_RETURN_VALUE; } // adjust min/max to match slab if (outExt[4] < slabMin) { outExt[4] = slabMin; } if (outExt[5] > slabMax) { outExt[5] = slabMax; } vtkCellLocator *locator = vtkCellLocator::New(); // Set up the cell locator. // If AutomaticOff, then NumberOfCellsPerBucket only used for allocating // memory. If AutomaticOn, then NumberOfCellsPerBucket is used to guess // the depth for the uniform octree required to support // NumberOfCellsPerBucket (assuming uniform distribution of cells). locator->SetDataSet( userData->Input[threadId] ); locator->AutomaticOff(); locator->SetMaxLevel( userData->Modeller->GetLocatorMaxLevel() ); locator->SetNumberOfCellsPerBucket( 1 ); locator->CacheCellBoundsOn(); locator->BuildLocator(); switch (userData->Modeller->GetOutputScalarType()) { vtkTemplateMacro( vtkMyImplicitModellerAppendExecute( userData->Modeller, userData->Input[threadId], output, outExt, userData->MaximumDistance, locator, threadId, static_cast(0))); default: vtkGenericWarningMacro("Execute: Unknown output ScalarType"); return VTK_THREAD_RETURN_VALUE; } locator->Delete(); return VTK_THREAD_RETURN_VALUE; } //---------------------------------------------------------------------------- // Templated append for VTK_CELL_MODE process mode and any type of output data template void vtkMyImplicitModellerAppendExecute(vtkMyImplicitModeller *self, vtkDataSet *input, vtkImageData *outData, double maxDistance, OT *) { int i, j, k, updateTime; vtkIdType cellNum; double *bounds, adjBounds[6]; double pcoords[3]; int outExt[6]; double x[3], prevDistance2, distance, distance2; int subId; double closestPoint[3]; double *weights=new double[input->GetMaxCellSize()]; double maxDistance2; double *spacing, *origin; spacing = outData->GetSpacing(); origin = outData->GetOrigin(); maxDistance2 = maxDistance * maxDistance; int *sampleDimensions = self->GetSampleDimensions(); // so we know how to scale if desired double scaleFactor = 0; // 0 used to indicate not scaling double toDoubleScaleFactor = 0; // 0 used to indicate not scaling double capValue = 0; // 0 used to indicate not clamping (float or double) if (self->GetOutputScalarType() != VTK_FLOAT && self->GetOutputScalarType() != VTK_DOUBLE) { capValue = self->GetCapValue(); if (self->GetScaleToMaximumDistance()) { scaleFactor = capValue / maxDistance; toDoubleScaleFactor = maxDistance / capValue; } } // // Traverse all cells; computing distance function on volume points. // vtkCell *cell; // MODIFICA: aggiunti sti 4 punti e il flag per calcolo dopo double p1[3], p2[3], p3[3], p4[3]; double tetravolume; int signFlag; updateTime = input->GetNumberOfCells() / 50; // update every 2% if (updateTime < 1) { updateTime = 1; } for (cellNum=0; cellNum < input->GetNumberOfCells(); cellNum++) { // MODIFICA: aggiunto signFlag = 1; cell = input->GetCell(cellNum); bounds = cell->GetBounds(); for (i=0; i<3; i++) { adjBounds[2*i] = bounds[2*i] - maxDistance; adjBounds[2*i+1] = bounds[2*i+1] + maxDistance; } // compute dimensional bounds in data set for (i = 0; i < 3; i++) { outExt[i*2] = (int) ((double)(adjBounds[2*i] - origin[i]) / spacing[i]); outExt[i*2 + 1] = (int) ((double)(adjBounds[2*i+1] - origin[i]) / spacing[i]); if (outExt[i*2] < 0) { outExt[i*2] = 0; } if (outExt[i*2 + 1] >= sampleDimensions[i]) { outExt[i*2 + 1] = sampleDimensions[i] - 1; } } vtkImageIterator outIt(outData, outExt); for (k = outExt[4]; k <= outExt[5]; k++) { x[2] = spacing[2] * k + origin[2]; for (j = outExt[2]; j <= outExt[3]; j++) { x[1] = spacing[1] * j + origin[1]; OT* outSI = outIt.BeginSpan(); for (i = outExt[0]; i <= outExt[1]; i++) { x[0] = spacing[0] * i + origin[0]; ConvertToDoubleDistance(*outSI, distance, prevDistance2, toDoubleScaleFactor); // union combination of distances // MODIFICARE QUI!!! if ( cell->EvaluatePosition(x, closestPoint, subId, pcoords, distance2, weights) != -1 && distance2 < prevDistance2 && distance2 <= maxDistance2 ) { // QUI SOTTO CALCOLARE SEGNO DEL VOLUME TETRAEDRO DEI PRIMI 3 PUNTI DELLA CELLA E DI X // E MOLTIPLICARVI LA DISTANCE, PRIMA CHE SIA PASSATA ALL OUTPUT // ATTENZIONE, FUNZIA SOLO NEL 3D E POTREBBE SBALLARE LE COSE SU MAXDISTANCE // USARE vtkMath::Determinant3x3 if ( cell->GetCellType() == VTK_TRIANGLE ) //dovrebbe render robusto ma imporre che il polydata passi per un trianglefilter { cell->GetPoints()->GetPoint(0,p4); cell->GetPoints()->GetPoint(1,p1); cell->GetPoints()->GetPoint(2,p2); tetravolume = vtkMath::Determinant3x3(p2[0]-p1[0], x[0]-p1[0], p4[0]-p1[0], p2[1]-p1[1], x[1]-p1[1], p4[1]-p1[1], p2[2]-p1[2], x[2]-p1[2], p4[2]-p1[2])/6.0; if (tetravolume < 0) signFlag = 1; else signFlag = -1; } // modificato col flag distance = sqrt(distance2)*signFlag; SetOutputDistance(distance, outSI, capValue, scaleFactor); } outSI++; } outIt.NextSpan(); } } if (cellNum % updateTime == 0) { self->UpdateProgress(double(cellNum + 1) / input->GetNumberOfCells()); } } delete [] weights; } // Append a data set to the existing output. To use this function, // you'll have to invoke the StartAppend() method before doing // successive appends. It's also a good idea to specify the model // bounds; otherwise the input model bounds is used. When you've // finished appending, use the EndAppend() method. void vtkMyImplicitModeller::Append(vtkDataSet *input) { vtkDebugMacro(<< "Appending data"); vtkImageData *output = this->GetOutput(); if ( !this->BoundsComputed ) { this->ComputeModelBounds(input); } if (this->ProcessMode == VTK_CELL_MODE) { if (!output->GetPointData()->GetScalars()) { vtkErrorMacro("Sanity check failed."); return; } switch (this->OutputScalarType) { vtkTemplateMacro( vtkMyImplicitModellerAppendExecute( this, input, output, this->InternalMaxDistance, static_cast(0))); } } else { vtkMyImplicitModellerAppendInfo info; double minZ, maxZ; int slabMin, slabMax, slabSize, i; vtkClipPolyData **minClipper = NULL, **maxClipper = NULL; vtkPlane ** minPlane = NULL, **maxPlane = NULL; double *spacing, *origin; spacing = output->GetSpacing(); origin = output->GetOrigin(); // Use a MultiThreader here, splitting the volume into slabs to be processed // by the separate threads // Set the number of threads to use, // then set the execution method and do it. this->Threader->SetNumberOfThreads( this->NumberOfThreads ); // set up the info object for the thread info.Modeller = this; info.MaximumDistance = this->InternalMaxDistance; info.Input = new vtkDataSet* [this->NumberOfThreads]; if (this->NumberOfThreads == 1) { info.Input[0] = input; } else { // if not PolyData, then copy the input for each thread if ( input->GetDataObjectType() != VTK_POLY_DATA ) { for (i = 0; i < this->NumberOfThreads; i++) { switch( input->GetDataObjectType() ) { case VTK_STRUCTURED_GRID: info.Input[i] = vtkStructuredGrid::New(); break; case VTK_IMAGE_DATA: info.Input[i] = vtkImageData::New(); break; case VTK_UNSTRUCTURED_GRID: info.Input[i] = vtkUnstructuredGrid::New(); break; case VTK_RECTILINEAR_GRID: info.Input[i] = vtkRectilinearGrid::New(); break; default: vtkErrorMacro(<<"Unexpected DataSet type!"); return; } info.Input[i]->CopyStructure(input); } } else // break up the input data into slabs to help ensure thread safety { minClipper = new vtkClipPolyData* [this->NumberOfThreads]; maxClipper = new vtkClipPolyData* [this->NumberOfThreads]; minPlane = new vtkPlane* [this->NumberOfThreads]; maxPlane = new vtkPlane* [this->NumberOfThreads]; slabSize = this->SampleDimensions[2] / this->NumberOfThreads; if (slabSize == 0) // in case threadCount > SampleDimensions[2] { slabSize = 1; } for (i = 0; i < this->NumberOfThreads; i++) { ////////////////////////////////////////////////// // do the 1st clip slabMin = i * slabSize; if (slabMin >= this->SampleDimensions[2]) { break; } // get/clip input cells in this slab + maxDistance+ minZ = spacing[2] * slabMin + origin[2] - this->InternalMaxDistance*1.00001; if (minZ < this->ModelBounds[4]) { minZ = this->ModelBounds[4]; } minPlane[i] = vtkPlane::New(); minPlane[i]->SetNormal(0.0, 0.0, -1.0); minPlane[i]->SetOrigin(0.0, 0.0, minZ); minClipper[i] = vtkClipPolyData::New(); minClipper[i]->SetInput((vtkPolyData *)input); minClipper[i]->SetClipFunction(minPlane[i]); minClipper[i]->SetValue( 0.0 ); minClipper[i]->InsideOutOn(); minClipper[i]->Update(); if ( minClipper[i]->GetOutput()->GetNumberOfCells() == 0 ) { info.Input[i] = NULL; maxPlane[i] = NULL; continue; } minClipper[i]->ReleaseDataFlagOn(); ////////////////////////////////////////////////// // do the 2nd clip slabMax = slabMin + slabSize - 1; if (i == this->NumberOfThreads - 1) { slabMax = this->SampleDimensions[2] - 1; } maxZ = spacing[2] * slabMax + origin[2] + this->InternalMaxDistance*1.00001; if (maxZ > this->ModelBounds[5]) { maxZ = this->ModelBounds[5]; } maxPlane[i] = vtkPlane::New(); maxPlane[i]->SetNormal(0.0, 0.0, 1.0); maxPlane[i]->SetOrigin(0.0, 0.0, maxZ); maxClipper[i] = vtkClipPolyData::New(); maxClipper[i]->SetInput(minClipper[i]->GetOutput()); maxClipper[i]->SetClipFunction(maxPlane[i]); maxClipper[i]->SetValue( 0.0 ); maxClipper[i]->InsideOutOn(); maxClipper[i]->Update(); if ( maxClipper[i]->GetOutput()->GetNumberOfCells() == 0 ) { info.Input[i] = NULL; } else { info.Input[i] = maxClipper[i]->GetOutput(); } } } } this->Threader->SetSingleMethod( vtkMyImplicitModeller_ThreadedAppend, (void *)&info); this->Threader->SingleMethodExecute(); // cleanup if (this->NumberOfThreads > 1) { if ( input->GetDataObjectType() != VTK_POLY_DATA ) { for (i = 0; i < this->NumberOfThreads; i++) { info.Input[i]->Delete(); } } else { for (i = 0; i < this->NumberOfThreads; i++) { minPlane[i]->Delete(); minClipper[i]->Delete(); if (maxPlane[i]) { maxPlane[i]->Delete(); maxClipper[i]->Delete(); } } delete [] minPlane; delete [] maxPlane; delete [] minClipper; delete [] maxClipper; } } delete [] info.Input; } } //---------------------------------------------------------------------------- // Method completes the append process (does the capping if requested). void vtkMyImplicitModeller::EndAppend() { vtkDataArray *newScalars; vtkDebugMacro(<< "End append"); if (!(newScalars =this->GetOutput()->GetPointData()->GetScalars())) { vtkErrorMacro("Sanity check failed."); return; } if ( this->Capping ) { this->Cap(newScalars); } this->UpdateProgress(1.0); } //---------------------------------------------------------------------------- int vtkMyImplicitModeller::RequestInformation ( vtkInformation * vtkNotUsed(request), vtkInformationVector ** vtkNotUsed( inputVector ), vtkInformationVector *outputVector) { // get the info objects vtkInformation* outInfo = outputVector->GetInformationObject(0); int i; double ar[3], origin[3]; vtkDataObject::SetPointDataActiveScalarInfo(outInfo, this->OutputScalarType, 1); outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), 0, this->SampleDimensions[0]-1, 0, this->SampleDimensions[1]-1, 0, this->SampleDimensions[2]-1); for (i=0; i < 3; i++) { origin[i] = this->ModelBounds[2*i]; if ( this->SampleDimensions[i] <= 1 ) { ar[i] = 1; } else { ar[i] = (this->ModelBounds[2*i+1] - this->ModelBounds[2*i]) / (this->SampleDimensions[i] - 1); } } outInfo->Set(vtkDataObject::ORIGIN(),origin,3); outInfo->Set(vtkDataObject::SPACING(),ar,3); return 1; } //---------------------------------------------------------------------------- int vtkMyImplicitModeller::RequestData( vtkInformation* vtkNotUsed( request ), vtkInformationVector** inputVector, vtkInformationVector* vtkNotUsed( outputVector )) { // get the input vtkInformation* inInfo = inputVector[0]->GetInformationObject(0); vtkDataSet *input = vtkDataSet::SafeDownCast( inInfo->Get(vtkDataObject::DATA_OBJECT())); vtkDebugMacro(<< "Executing implicit model"); if (input == NULL) { // we do not want to release the data because user might // have called Append ... return 0; } this->StartAppend(1); this->Append(input); this->EndAppend(); return 1; } // Compute ModelBounds from input geometry. double vtkMyImplicitModeller::ComputeModelBounds(vtkDataSet *input) { double *bounds, maxDist; int i; vtkImageData *output=this->GetOutput(); double tempd[3]; // compute model bounds if not set previously if ( this->ModelBounds[0] >= this->ModelBounds[1] || this->ModelBounds[2] >= this->ModelBounds[3] || this->ModelBounds[4] >= this->ModelBounds[5] ) { if (input != NULL) { bounds = input->GetBounds(); } else { vtkDataSet *dsInput = vtkDataSet::SafeDownCast(this->GetInput()); if (dsInput != NULL) { bounds = dsInput->GetBounds(); } else { vtkErrorMacro( << "An input must be specified to Compute the model bounds."); return VTK_FLOAT_MAX; } } } else { bounds = this->ModelBounds; } for (maxDist=0.0, i=0; i<3; i++) { if ( (bounds[2*i+1] - bounds[2*i]) > maxDist ) { maxDist = bounds[2*i+1] - bounds[2*i]; } } // adjust bounds so model fits strictly inside (only if not set previously) if ( this->AdjustBounds ) { for (i=0; i<3; i++) { this->ModelBounds[2*i] = bounds[2*i] - maxDist*this->AdjustDistance; this->ModelBounds[2*i+1] = bounds[2*i+1] + maxDist*this->AdjustDistance; } } else // to handle problem case where bounds not specified and AdjustBounds // not on; will be setting ModelBounds to self if previosusly set { for (i=0; i<3; i++) { this->ModelBounds[2*i] = bounds[2*i]; this->ModelBounds[2*i+1] = bounds[2*i+1]; } } maxDist *= this->MaximumDistance; // Set volume origin and data spacing output->SetOrigin(this->ModelBounds[0], this->ModelBounds[2], this->ModelBounds[4]); for (i=0; i<3; i++) { tempd[i] = (this->ModelBounds[2*i+1] - this->ModelBounds[2*i]) / (this->SampleDimensions[i] - 1); } output->SetSpacing(tempd); vtkInformation *outInfo = this->GetExecutive()->GetOutputInformation(0); outInfo->Set(vtkDataObject::ORIGIN(),this->ModelBounds[0], this->ModelBounds[2], this->ModelBounds[4]); outInfo->Set(vtkDataObject::SPACING(),tempd,3); this->BoundsComputed = 1; this->InternalMaxDistance = maxDist; return maxDist; } //---------------------------------------------------------------------------- // Set the i-j-k dimensions on which to sample the distance function. void vtkMyImplicitModeller::SetSampleDimensions(int i, int j, int k) { int dim[3]; dim[0] = i; dim[1] = j; dim[2] = k; this->SetSampleDimensions(dim); } //---------------------------------------------------------------------------- void vtkMyImplicitModeller::SetSampleDimensions(int dim[3]) { int dataDim, i; vtkDebugMacro(<< " setting SampleDimensions to (" << dim[0] << "," << dim[1] << "," << dim[2] << ")"); if ( dim[0] != this->SampleDimensions[0] || dim[1] != this->SampleDimensions[1] || dim[2] != this->SampleDimensions[2] ) { if ( dim[0]<1 || dim[1]<1 || dim[2]<1 ) { vtkErrorMacro (<< "Bad Sample Dimensions, retaining previous values"); return; } for (dataDim=0, i=0; i<3 ; i++) { if (dim[i] > 1) { dataDim++; } } if ( dataDim < 3 ) { vtkErrorMacro(<<"Sample dimensions must define a volume!"); return; } for ( i=0; i<3; i++) { this->SampleDimensions[i] = dim[i]; } this->Modified(); } } //---------------------------------------------------------------------------- void vtkMyImplicitModeller::Cap(vtkDataArray *s) { int i,j,k; int idx; int d01=this->SampleDimensions[0]*this->SampleDimensions[1]; // i-j planes k = 0; for (j=0; jSampleDimensions[1]; j++) { for (i=0; iSampleDimensions[0]; i++) { s->SetComponent(i+j*this->SampleDimensions[0],0, this->CapValue); } } k = this->SampleDimensions[2] - 1; idx = k*d01; for (j=0; jSampleDimensions[1]; j++) { for (i=0; iSampleDimensions[0]; i++) { s->SetComponent(idx+i+j*this->SampleDimensions[0], 0, this->CapValue); } } // j-k planes i = 0; for (k=0; kSampleDimensions[2]; k++) { for (j=0; jSampleDimensions[1]; j++) { s->SetComponent(j*this->SampleDimensions[0]+k*d01,0,this->CapValue); } } i = this->SampleDimensions[0] - 1; for (k=0; kSampleDimensions[2]; k++) { for (j=0; jSampleDimensions[1]; j++) { s->SetComponent(i+j*this->SampleDimensions[0]+k*d01,0, this->CapValue); } } // i-k planes j = 0; for (k=0; kSampleDimensions[2]; k++) { for (i=0; iSampleDimensions[0]; i++) { s->SetComponent(i+k*d01,0, this->CapValue); } } j = this->SampleDimensions[1] - 1; idx = j*this->SampleDimensions[0]; for (k=0; kSampleDimensions[2]; k++) { for (i=0; iSampleDimensions[0]; i++) { s->SetComponent(idx+i+k*d01,0, this->CapValue); } } } //---------------------------------------------------------------------------- const char *vtkMyImplicitModeller::GetProcessModeAsString() { if (this->ProcessMode == VTK_CELL_MODE) { return "PerCell"; } else { return "PerVoxel"; } } //---------------------------------------------------------------------------- int vtkMyImplicitModeller::FillInputPortInformation( int vtkNotUsed( port ), vtkInformation* info) { info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataSet"); info->Set(vtkAlgorithm::INPUT_IS_OPTIONAL(), 1); return 1; } //---------------------------------------------------------------------------- int vtkMyImplicitModeller::ProcessRequest(vtkInformation* request, vtkInformationVector** inputVector, vtkInformationVector* outputVector) { // If we have no input then we will not generate the output because // the user already called StartAppend/Append/EndAppend. if(request->Has(vtkDemandDrivenPipeline::REQUEST_DATA_NOT_GENERATED())) { if(inputVector[0]->GetNumberOfInformationObjects() == 0) { vtkInformation* outInfo = outputVector->GetInformationObject(0); outInfo->Set(vtkDemandDrivenPipeline::DATA_NOT_GENERATED(), 1); } return 1; } else if(request->Has(vtkDemandDrivenPipeline::REQUEST_DATA())) { if(inputVector[0]->GetNumberOfInformationObjects() == 0) { return 1; } } return this->Superclass::ProcessRequest(request, inputVector, outputVector); } void vtkMyImplicitModeller::PrintSelf(ostream& os, vtkIndent indent) { this->Superclass::PrintSelf(os,indent); os << indent << "Maximum Distance: " << this->MaximumDistance << "\n"; os << indent << "OutputScalarType: " << this->OutputScalarType << "\n"; os << indent << "Sample Dimensions: (" << this->SampleDimensions[0] << ", " << this->SampleDimensions[1] << ", " << this->SampleDimensions[2] << ")\n"; os << indent << "ModelBounds: \n"; os << indent << " Xmin,Xmax: (" << this->ModelBounds[0] << ", " << this->ModelBounds[1] << ")\n"; os << indent << " Ymin,Ymax: (" << this->ModelBounds[2] << ", " << this->ModelBounds[3] << ")\n"; os << indent << " Zmin,Zmax: (" << this->ModelBounds[4] << ", " << this->ModelBounds[5] << ")\n"; os << indent << "ScaleToMaximumDistance: " << (this->ScaleToMaximumDistance ? "On\n" : "Off\n"); os << indent << "AdjustBounds: " << (this->AdjustBounds ? "On\n" : "Off\n"); os << indent << "Adjust Distance: " << this->AdjustDistance << "\n"; os << indent << "Process Mode: " << this->ProcessMode << "\n"; os << indent << "Locator Max Level: " << this->LocatorMaxLevel << "\n"; os << indent << "Capping: " << (this->Capping ? "On\n" : "Off\n"); os << indent << "Cap Value: " << this->CapValue << "\n"; os << indent << "Process Mode: " << this->GetProcessModeAsString() << endl; os << indent << "Number Of Threads (for PerVoxel mode): " << this->NumberOfThreads << endl; }