41 #include <vtkImageData.h>
61 std::vector<PropertyPtr> retval;
62 retval.push_back(this->getInterpolationStepsOption(root));
66 DoublePropertyPtr PNNReconstructionMethodService::getInterpolationStepsOption(QDomElement root)
70 "Interpolation steps in voxels", 3,
DoubleRange(1, 10, 1), 0, root);
76 double* t = tt.begin();
80 (*p)[0] = t[0] * x + t[1] * y + t[2] * z + t[3];
81 (*p)[1] = t[4] * x + t[5] * y + t[6] * z + t[7];
82 (*p)[2] = t[8] * x + t[9] * y + t[10] * z + t[11];
90 std::vector<TimedPosition> frameInfo = input->getFrames();
91 if (frameInfo.empty())
93 if (input->getDimensions()[2]==0)
98 Eigen::Array3i inputDims = input->getDimensions();
100 Eigen::Array3i targetDims(target->GetDimensions());
101 Vector3D targetSpacing(target->GetSpacing());
107 int* outputDims = tempOutput->GetDimensions();
110 if (inputDims[2] != static_cast<int> (frameInfo.size()))
114 Vector3D inputSpacing(input->getSpacing());
115 Vector3D outputSpacing(tempOutput->GetSpacing());
118 unsigned char *outputPointer =
static_cast<unsigned char*
> (tempOutput->GetScalarPointer());
119 unsigned char* maskPointer =
static_cast<unsigned char*
> (input->getMask()->GetScalarPointer());
122 for (
int record = 0; record < inputDims[2]; record++)
124 unsigned char *inputPointer = input->getFrame(record);
125 boost::array<double, 16> recordTransform = frameInfo[record].mPos.flatten();
127 for (
int beam = 0; beam < inputDims[0]; beam++)
129 for (
int sample = 0; sample < inputDims[1]; sample++)
131 if (!validPixel(beam, sample, inputDims, maskPointer))
133 Vector3D inputPoint(beam * inputSpacing[0], sample * inputSpacing[1], 0.0);
136 int outputVoxelX =
static_cast<int> ((outputPoint[0] / outputSpacing[0]) + 0.5);
137 int outputVoxelY =
static_cast<int> ((outputPoint[1] / outputSpacing[1]) + 0.5);
138 int outputVoxelZ =
static_cast<int> ((outputPoint[2] / outputSpacing[2]) + 0.5);
140 if (validVoxel(outputVoxelX, outputVoxelY, outputVoxelZ, outputDims))
142 int outputIndex = outputVoxelX + outputVoxelY * outputDims[0] + outputVoxelZ * outputDims[0]
144 int inputIndex = beam + sample * inputDims[0];
147 outputPointer[outputIndex] = std::max<unsigned char>(inputPointer[inputIndex], outputPointer[outputIndex]);
149 outputPointer[outputIndex] = std::max<unsigned char>(inputPointer[inputIndex], 1);
157 this->interpolate(tempOutputData, outputData, settings);
167 inline int getIndex_z_last(
int x,
int y,
int z,
const Eigen::Array3i& dim)
169 return x + y*dim[0] + z*dim[0]*dim[1];
173 inline int getIndex_x_last(
int y,
int z,
int x,
const Eigen::Array3i& dim)
175 return x + y*dim[0] + z*dim[0]*dim[1];
179 inline int getIndex_y_last(
int z,
int x,
int y,
const Eigen::Array3i& dim)
181 return x + y*dim[0] + z*dim[0]*dim[1];
190 template <
class FUNCTION>
191 void maskAlongDim(
int a_dim,
int b_dim,
int c_dim,
const Eigen::Array3i& dim,
unsigned char *inputPtr,
unsigned char *maskPtr, FUNCTION getIndex)
193 for (
int a = 0; a < a_dim; a++)
195 for (
int b = 0; b < b_dim; b++)
199 for (
int c = 0; c < c_dim; c++)
201 int index = getIndex(a, b, c, dim);
202 if (inputPtr[index]>0)
208 for (
int c = c_dim-1; c >=0; c--)
210 int index = getIndex(a, b, c, dim);
211 if (inputPtr[index]>0)
217 for (
int c = start; c <= stop; c++)
219 int index = getIndex(a, b, c, dim);
237 Eigen::Array3i dim(inputData->GetDimensions());
238 Vector3D spacing(inputData->GetSpacing());
240 unsigned char *inputPtr =
static_cast<unsigned char*
> (inputData->GetScalarPointer());
241 unsigned char *maskPtr =
static_cast<unsigned char*
> (mask->GetScalarPointer());
244 maskAlongDim(dim[0], dim[1], dim[2], dim, inputPtr, maskPtr, &getIndex_z_last);
245 maskAlongDim(dim[1], dim[2], dim[0], dim, inputPtr, maskPtr, &getIndex_x_last);
246 maskAlongDim(dim[2], dim[0], dim[1], dim, inputPtr, maskPtr, &getIndex_y_last);
254 DoublePropertyPtr interpolationStepsOption = this->getInterpolationStepsOption(settings);
255 int interpolationSteps =
static_cast<int> (interpolationStepsOption->getValue());
261 Eigen::Array3i outputDims(output->GetDimensions());
263 Eigen::Array3i inputDims(input->GetDimensions());
265 unsigned char *inputPointer =
static_cast<unsigned char*
> (input->GetScalarPointer());
266 unsigned char *outputPointer =
static_cast<unsigned char*
> (output->GetScalarPointer());
267 unsigned char *maskPointer =
static_cast<unsigned char*
> (mask->GetScalarPointer());
269 if ((outputDims[0] != inputDims[0]) || (outputDims[1] != inputDims[1]) || (outputDims[2] != inputDims[2]))
274 int total = outputDims[0] * outputDims[1] * outputDims[2];
278 for (
int x = 0; x < outputDims[0]; x++)
280 for (
int y = 0; y < outputDims[1]; y++)
282 for (
int z = 0; z < outputDims[2]; z++)
284 int outputIndex = x + y * outputDims[0] + z * outputDims[0] * outputDims[1];
287 if (maskPointer[outputIndex]==0)
292 else if (inputPointer[outputIndex]>0)
294 outputPointer[outputIndex] = inputPointer[outputIndex];
300 this->fillHole(inputPointer, outputPointer, x, y, z, outputDims, interpolationSteps);
306 int valid = 100*double(ignored)/double(total);
307 int outside = 100*double(removed)/double(total);
308 int holes = 100*double(total-ignored-removed)/double(total);
310 QString(
"PNN: Size: %1Mb, Valid voxels: %2\%, Outside mask: %3\% Filled holes [steps=%4, %5s]: %6\%")
311 .arg(total/1024/1024)
314 .arg(interpolationSteps)
315 .arg(timer.getElapsedSecondsAsString())
323 void PNNReconstructionMethodService::fillHole(
unsigned char *inputPointer,
unsigned char *outputPointer,
int x,
int y,
int z,
const Eigen::Array3i& dim,
int interpolationSteps)
325 int outputIndex = x + y * dim[0] + z * dim[0] * dim[1];
326 bool interpolated =
false;
334 for (
int i = -localArea; i < localArea + 1; i++)
336 for (
int j = -localArea; j < localArea + 1; j++)
338 for (
int k = -localArea; k < localArea + 1; k++)
340 int localIndex = outputIndex + i + j*dim[0] + k*dim[0]*dim[1];
342 if (validVoxel(x + i, y + j, z + k, dim.data()) && inputPointer[localIndex] > 0.1)
344 tempVal += inputPointer[localIndex];
360 outputPointer[outputIndex] =
static_cast<int> ((tempVal / count) + 0.5);
361 outputPointer[outputIndex] = std::max<unsigned char>(1, outputPointer[outputIndex]);
367 }
while (localArea <= interpolationSteps && !interpolated);
QString qstring_cast(const T &val)
Utility class for describing a bounded numeric range.
boost::shared_ptr< class Image > ImagePtr
virtual QString getName() const
virtual bool reconstruct(ProcessedUSInputDataPtr input, vtkImageDataPtr outputData, QDomElement settings)
void reportWarning(QString msg)
Settings * settings()
Shortcut for accessing the settings instance.
vtkImageDataPtr generateVtkImageData(Eigen::Array3i dim, Vector3D spacing, const unsigned char initValue, int components)
PNNReconstructionMethodService(ctkPluginContext *context)
virtual std::vector< PropertyPtr > getSettings(QDomElement root)
virtual ~PNNReconstructionMethodService()
boost::shared_ptr< class DoubleProperty > DoublePropertyPtr
Eigen::Vector3d Vector3D
Vector3D is a representation of a point or vector in 3D.
void setDeepModified(vtkImageDataPtr image)
static DoublePropertyPtr initialize(const QString &uid, QString name, QString help, double value, DoubleRange range, int decimals, QDomNode root=QDomNode())
boost::shared_ptr< class ProcessedUSInputData > ProcessedUSInputDataPtr
void optimizedCoordTransform(Vector3D *p, boost::array< double, 16 > tt)
void reportDebug(QString msg)