<commit_before>// --------------------------------------------------------------------------
// File:             SparseToDenseFieldResampleFilter.hxx
// Date:             30/10/2012
// Author:           <EMAIL> (<NAME>, OE)
// Version:          0.1
// License:          BSD
// --------------------------------------------------------------------------
//
// Copyright (c) 2012, <EMAIL> (<NAME>)
// with Signal Processing Lab 5, EPFL (LTS5-EPFL)
// and Biomedical Image Technology, UPM (BIT-UPM)
// All rights reserved.
// 
// This file is part of ACWEReg
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
// * Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
// * Neither the names of the LTS5-EFPL and the BIT-UPM, nor the names of its
// contributors may be used to endorse or promote products derived from this
// software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY Oscar Esteban ''AS IS'' AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
// DISCLAIMED. IN NO EVENT SHALL OSCAR ESTEBAN BE LIABLE FOR ANY
// DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
// LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
// ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

#ifndef SPARSETODENSEFIELDRESAMPLEFILTER_HXX_
#define SPARSETODENSEFIELDRESAMPLEFILTER_HXX_

#include "SparseToDenseFieldResampleFilter.h"
#include "IDWMultivariateInterpolator.h"
#include "itkImageRegionIteratorWithIndex.h"

namespace rstk {

template<class TInputMesh, class TOutputImage>
SparseToDenseFieldResampleFilter<TInputMesh, TOutputImage>::SparseToDenseFieldResampleFilter() {
	m_OutputSpacing.Fill(1.0);
	m_OutputOrigin.Fill(0.0);
	m_OutputDirection.SetIdentity();
	m_OutputSize.Fill(0);
	m_OutputStartIndex.Fill(0);

	m_Interpolator =
			IDWMultivariateInterpolator<InputMeshType, OutputImageType>::New();
	m_DefaultPixelValue = NumericTraits<ValueType>::ZeroValue(
			m_DefaultPixelValue);
	m_IsPhiInitialized = false;
	m_N = 0;
	m_k = 0;
};

template<class TInputMesh, class TOutputImage>
void SparseToDenseFieldResampleFilter<TInputMesh, TOutputImage>::
AddControlPoints( const typename SparseToDenseFieldResampleFilter<TInputMesh, TOutputImage>::InputMeshType* prior) {
	for( size_t i = 0; i<prior->GetNumberOfPoints();i++ ) {
		this->m_ControlPoints.push_back( prior->GetPoint(i));
	}
	this->m_N = this->m_ControlPoints.size();
}

template<class TInputMesh, class TOutputImage>
void SparseToDenseFieldResampleFilter<TInputMesh, TOutputImage>::PrintSelf(
		std::ostream &os, Indent indent) const {
	//Superclass::PrintSelf(os, indent);

	os << indent << "DefaultPixelValue: "
			<< static_cast<typename NumericTraits<OutputPixelType>::PrintType>(m_DefaultPixelValue)
			<< std::endl;
	os << indent << "Size: " << m_OutputSize << std::endl;
	os << indent << "OutputStartIndex: " << m_OutputStartIndex << std::endl;
	os << indent << "OutputSpacing: " << m_OutputSpacing << std::endl;
	os << indent << "OutputOrigin: " << m_OutputOrigin << std::endl;
	os << indent << "OutputDirection: " << m_OutputDirection << std::endl;
	os << indent << "Interpolator: " << m_Interpolator.GetPointer()
			<< std::endl;
}
;

/**
 * Set the output image spacing.
 */
template<class TInputMesh, class TOutputImage>
void SparseToDenseFieldResampleFilter<TInputMesh, TOutputImage>::SetOutputSpacing(
		const double *spacing) {
	OutputSpacingType s(spacing);

	this->SetOutputSpacing(s);
}

/**
 * Set the output image origin.
 */
template<class TInputMesh, class TOutputImage>
void SparseToDenseFieldResampleFilter<TInputMesh, TOutputImage>::SetOutputOrigin(
		const double *origin) {
	OutputPointType p(origin);
	this->SetOutputOrigin(p);
}

/**
 * Set up state of filter before multi-threading.
 * InterpolatorType::SetInputImage is not thread-safe and hence
 * has to be set up before ThreadedGenerateData
 */
template<class TInputMesh, class TOutputImage>
void SparseToDenseFieldResampleFilter<TInputMesh, TOutputImage>::GenerateData() {
	if ( m_N == 0 ) {
		itkExceptionMacro(<< "Shape priors not set");
	}

	if (m_DefaultPixelValue.Size() == 0) {
		NumericTraits<OutputPixelType>::SetLength(m_DefaultPixelValue,
				this->GetOutput()->GetNumberOfComponentsPerPixel());
		m_DefaultPixelValue.Fill(0);
	}

//	this->GetOutput()->FillBuffer( this->m_DefaultPixelValue );

	m_k = this->GetOutput()->GetLargestPossibleRegion().GetNumberOfPixels();

	for (size_t i = 0; i < OutputImageDimension; i++ ) m_LevelSetVector[i] = SpeedsVector( m_N );
	size_t nConts = this->GetNumberOfIndexedInputs();

	if (!m_IsPhiInitialized) {
		for( size_t cont = 0; cont< nConts; cont++ ) {
			m_Phi = WeightsMatrix(m_k, m_N);
			// Walk the output region
			size_t row, col;
			double dist, wi;
			OutputPointType ci, pi;

			for( row = 0; row < this->m_k; row++ ) {
				this->GetOutput()->TransformIndexToPhysicalPoint(this->GetOutput()->ComputeIndex(row), pi);

				for (col = 0; col < this->m_N; col++) {
					// TODO Extract RBF from here (use a functor?)
					// Compute weight i: wi(x) (inverse distance)
					dist = (this->m_ControlPoints[col] - pi).GetNorm();
					wi = (dist < vnl_math::eps) ? 1.0 : (1.0 / vcl_pow(dist, 2));
					if (wi > 1e-3) {
						m_Phi.put(row, col, wi);
					}
				}
			}
		}
		m_IsPhiInitialized = true;
	}

	size_t contsOffset = 0;
	for( size_t cont = 0; cont< nConts; cont++ ) {
		OutputPixelType ni;
		size_t nPoints = this->GetInput(cont)->GetNumberOfPoints();
		for (size_t el = 0; el < nPoints; el++) {
			this->GetInput(cont)->GetPointData( el, &ni );

			for ( size_t i = 0; i < OutputImageDimension; i++ ) {
				if( std::isnan( ni[i] )) ni[i] = 0;
				m_LevelSetVector[i].put( contsOffset + el, ni[i] );
			}
		}
		contsOffset += nPoints;
	}

	double norms = 0.0;
	for (size_t i = 0; i < OutputImageDimension; i++ ) {
		norms+= m_LevelSetVector[i].one_norm();
	}


	if( norms==0 ) return;

    SpeedsVector outvector[3];
    norms=0;
    for ( size_t i = 0; i < OutputImageDimension; i++ ) {
    	m_Phi.mult(m_LevelSetVector[i], outvector[i] );
    	norms+= outvector[i].one_norm();<commit_msg>Fixed #45: contours iterated twice<commit_after>// --------------------------------------------------------------------------
// File:             SparseToDenseFieldResampleFilter.hxx
// Date:             30/10/2012
// Author:           <EMAIL> (<NAME>, OE)
// Version:          0.1
// License:          BSD
// --------------------------------------------------------------------------
//
// Copyright (c) 2012, <EMAIL> (<NAME>)
// with Signal Processing Lab 5, EPFL (LTS5-EPFL)
// and Biomedical Image Technology, UPM (BIT-UPM)
// All rights reserved.
// 
// This file is part of ACWEReg
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
// * Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
// * Neither the names of the LTS5-EFPL and the BIT-UPM, nor the names of its
// contributors may be used to endorse or promote products derived from this
// software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY Oscar Esteban ''AS IS'' AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
// DISCLAIMED. IN NO EVENT SHALL OSCAR ESTEBAN BE LIABLE FOR ANY
// DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
// LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
// ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

#ifndef SPARSETODENSEFIELDRESAMPLEFILTER_HXX_
#define SPARSETODENSEFIELDRESAMPLEFILTER_HXX_

#include "SparseToDenseFieldResampleFilter.h"
#include "IDWMultivariateInterpolator.h"
#include "itkImageRegionIteratorWithIndex.h"

namespace rstk {

template<class TInputMesh, class TOutputImage>
SparseToDenseFieldResampleFilter<TInputMesh, TOutputImage>::SparseToDenseFieldResampleFilter() {
	m_OutputSpacing.Fill(1.0);
	m_OutputOrigin.Fill(0.0);
	m_OutputDirection.SetIdentity();
	m_OutputSize.Fill(0);
	m_OutputStartIndex.Fill(0);

	m_Interpolator =
			IDWMultivariateInterpolator<InputMeshType, OutputImageType>::New();
	m_DefaultPixelValue = NumericTraits<ValueType>::ZeroValue(
			m_DefaultPixelValue);
	m_IsPhiInitialized = false;
	m_N = 0;
	m_k = 0;
};

template<class TInputMesh, class TOutputImage>
void SparseToDenseFieldResampleFilter<TInputMesh, TOutputImage>::
AddControlPoints( const typename SparseToDenseFieldResampleFilter<TInputMesh, TOutputImage>::InputMeshType* prior) {
	for( size_t i = 0; i<prior->GetNumberOfPoints();i++ ) {
		this->m_ControlPoints.push_back( prior->GetPoint(i));
	}
	this->m_N = this->m_ControlPoints.size();
}

template<class TInputMesh, class TOutputImage>
void SparseToDenseFieldResampleFilter<TInputMesh, TOutputImage>::PrintSelf(
		std::ostream &os, Indent indent) const {
	//Superclass::PrintSelf(os, indent);

	os << indent << "DefaultPixelValue: "
			<< static_cast<typename NumericTraits<OutputPixelType>::PrintType>(m_DefaultPixelValue)
			<< std::endl;
	os << indent << "Size: " << m_OutputSize << std::endl;
	os << indent << "OutputStartIndex: " << m_OutputStartIndex << std::endl;
	os << indent << "OutputSpacing: " << m_OutputSpacing << std::endl;
	os << indent << "OutputOrigin: " << m_OutputOrigin << std::endl;
	os << indent << "OutputDirection: " << m_OutputDirection << std::endl;
	os << indent << "Interpolator: " << m_Interpolator.GetPointer()
			<< std::endl;
}
;

/**
 * Set the output image spacing.
 */
template<class TInputMesh, class TOutputImage>
void SparseToDenseFieldResampleFilter<TInputMesh, TOutputImage>::SetOutputSpacing(
		const double *spacing) {
	OutputSpacingType s(spacing);

	this->SetOutputSpacing(s);
}

/**
 * Set the output image origin.
 */
template<class TInputMesh, class TOutputImage>
void SparseToDenseFieldResampleFilter<TInputMesh, TOutputImage>::SetOutputOrigin(
		const double *origin) {
	OutputPointType p(origin);
	this->SetOutputOrigin(p);
}

/**
 * Set up state of filter before multi-threading.
 * InterpolatorType::SetInputImage is not thread-safe and hence
 * has to be set up before ThreadedGenerateData
 */
template<class TInputMesh, class TOutputImage>
void SparseToDenseFieldResampleFilter<TInputMesh, TOutputImage>::GenerateData() {
	if ( m_N == 0 ) {
		itkExceptionMacro(<< "Shape priors not set");
	}

	if (m_DefaultPixelValue.Size() == 0) {
		NumericTraits<OutputPixelType>::SetLength(m_DefaultPixelValue,
				this->GetOutput()->GetNumberOfComponentsPerPixel());
		m_DefaultPixelValue.Fill(0);
	}

//	this->GetOutput()->FillBuffer( this->m_DefaultPixelValue );

	m_k = this->GetOutput()->GetLargestPossibleRegion().GetNumberOfPixels();

	for (size_t i = 0; i < OutputImageDimension; i++ ) m_LevelSetVector[i] = SpeedsVector( m_N );
	size_t nConts = this->GetNumberOfIndexedInputs();

	if (!m_IsPhiInitialized) {
		m_Phi = WeightsMatrix(m_k, m_N);
		double dist, wi;
		OutputPointType ci, pi;
		size_t row, col;

		// Walk the output region
		for( row = 0; row < this->m_k; row++ ) {
			this->GetOutput()->TransformIndexToPhysicalPoint(this->GetOutput()->ComputeIndex(row), pi);

			for ( col=0; col < this->m_N; col++) {
				// TODO Extract RBF from here (use a functor?)
				// Compute weight i: wi(x) (inverse distance)
				dist = (this->m_ControlPoints[col] - pi).GetNorm();
				wi = (dist < vnl_math::eps) ? 1.0 : (1.0 / vcl_pow(dist, 2));
				if (wi > 1e-3) {
					m_Phi.put(row, col, wi);
				}
			}
		}

		m_IsPhiInitialized = true;
	}

	size_t contsOffset = 0;
	for( size_t cont = 0; cont< nConts; cont++ ) {
		OutputPixelType ni;
		size_t nPoints = this->GetInput(cont)->GetNumberOfPoints();
		for (size_t el = 0; el < nPoints; el++) {
			this->GetInput(cont)->GetPointData( el, &ni );

			for ( size_t i = 0; i < OutputImageDimension; i++ ) {
				if( std::isnan( ni[i] )) ni[i] = 0;
				m_LevelSetVector[i].put( contsOffset + el, ni[i] );
			}
		}
		contsOffset += nPoints;
	}

	double norms = 0.0;
	for (size_t i = 0; i < OutputImageDimension; i++ ) {
		norms+= m_LevelSetVector[i].one_norm();
	}


	if( norms==0 ) return;

    SpeedsVector outvector[3];
    norms=0;
    for ( size_t i = 0; i < OutputImageDimension; i++ ) {
    	m_Phi.mult(m_LevelSetVector[i], outvector[i] );
    	norms+= outvector[i].one_norm();
{'fraction_non_alphanumeric': 0.14454819498566282, 'fraction_numerical': 0.007720020586721565, 'mean_word_length': 3.9805931893079456, 'pattern_counts': {'":': 0, '<': 184, '<?xml version=': 0, '>': 101, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '979928', 'n_tokens_mistral': 4676, 'n_tokens_neox': 4305, 'n_words': 1181}