<commit_before>/// @file mist/registration.h
//!
//! @brief ́E񍄑̃WXg[V
//!
#ifndef __INCLUDE_MIST_REGISTRATION__
#define __INCLUDE_MIST_REGISTRATION__


#ifndef __INCLUDE_MIST_CONF_H__
#include "config/mist_conf.h"
#endif


#ifndef __INCLUDE_MIST_H__
#include "mist.h"
#endif

#ifndef __INCLUDE_MIST_MINIMIZATION__
#include "minimization.h"
#endif

#ifndef __INCLUDE_MIST_DRAWING__
#include "drawing.h"
#endif


// mistOԂ̎n܂
_MIST_BEGIN

namespace __non_rigid_registration_utility__
{
	inline void FFD( double v, double &f0, double &f1, double &f2, double &f3 )
	{
		double v_ = 1.0 - v;
		double v2 = v * v;
		double v3 = v2 * v;
		f0 = v_ * v_ * v_ / 6.0;
		f1 = 0.5 * v3 - v2 + 4.0 / 6.0;
		f2 = ( v2 - v3 + v ) * 0.5 + 1.0 / 6.0;
		f3 = v3 / 6.0;
	}


	// _pāC\[X摜^[Qbg摜ɕό`
	template < class TARGETTYPE, class TARGETTYPEA, class SOURCETYPE, class SOURCETYPEA, class CONTROLMESH >
	static void non_rigid_transformation_( array2< TARGETTYPE, TARGETTYPEA > &target, const array2< SOURCETYPE, SOURCETYPEA > &source, const CONTROLMESH &control_mesh,
											typename CONTROLMESH::difference_type mx = -1,
											typename CONTROLMESH::difference_type my = -1,
											typename CONTROLMESH::difference_type mz = -1,
											typename CONTROLMESH::size_type thread_id = 0,
											typename CONTROLMESH::size_type thread_num = 1
										)
	{
		typedef array2< TARGETTYPE, TARGETTYPEA >  target_image_type;
		typedef array2< SOURCETYPE, SOURCETYPEA >  source_image_type;
		typedef CONTROLMESH control_mesh_type;

		typedef typename target_image_type::size_type size_type;
		typedef typename target_image_type::difference_type difference_type;

		difference_type i, j;
		double u, v;
		double sx[ 4 ], sy[ 4 ];
		double stepW = control_mesh.width( ) == 1 ? 1.0 : target.width( ) / static_cast< double >( control_mesh.width( ) - 1 );
		double stepH = control_mesh.height( ) == 1 ? 1.0 : target.height( ) / static_cast< double >( control_mesh.height( ) - 1 );
		double _1_stepW = 1.0 / stepW;
		double _1_stepH = 1.0 / stepH;
		double _1_ax = 1.0 / source.reso1( );
		double _1_ay = 1.0 / source.reso2( );

		//stepW *= target.reso1( );
		//stepH *= target.reso2( );

		difference_type d0, d1, d2, d3;
		{
			difference_type cx = source.width( ) / 2;
			difference_type cy = source.height( ) / 2;
			typename source_image_type::const_pointer ppp = &source( cx, cy );
			d0 = 0;
			d1 = &source( cx    , cy + 1 ) - ppp;
			d2 = &source( cx + 1, cy + 1 ) - ppp;
			d3 = &source( cx + 1, cy     ) - ppp;
		}

		difference_type isx, isy, iex, iey;
		difference_type mw = control_mesh.width( );
		difference_type mh = control_mesh.height( );
		if( mx < 0 || mw <= mx || my < 0 || mh <= my )
		{
			isx = 0;
			iex = target.width( ) - 1;
			isy = 0;
			iey = target.height( ) - 1;
		}
		else
		{
			isx = static_cast< difference_type >( ( mx - 2 ) * stepW );
			iex = static_cast< difference_type >( ( mx + 2 ) * stepW );
			isy = static_cast< difference_type >( ( my - 2 ) * stepH );
			iey = static_cast< difference_type >( ( my + 2 ) * stepH );

			difference_type tw = target.width( );
			difference_type th = target.height( );
			isx = isx > 0  ? isx : 0;
			isy = isy > 0  ? isy : 0;
			iex = iex < tw ? iex : tw - 1;
			iey = iey < th ? iey : th - 1;
		}


		for( difference_type y = isy + thread_id ; y <= iey ; y += thread_num )
		{
			v = y * _1_stepH;
			j = static_cast< difference_type >( v );
			v -= j;
			j--;

			FFD( v, sy[ 0 ], sy[ 1 ], sy[ 2 ], sy[ 3 ] );

			for( difference_type x = isx ; x <= iex ; x++ )
			{
				u = x * _1_stepW;
				i = static_cast< difference_type >( u );
				u -= i;
				i--;

				FFD( u, sx[ 0 ], sx[ 1 ], sx[ 2 ], sx[ 3 ] );

				double xx = 0.0, yy = 0.0;
				for( difference_type m = 0 ; m <= 3 ; m++ )
				{
					typename control_mesh_type::const_pointer p = &control_mesh( i, j + m, 0 );
					for( difference_type l = 0 ; l <= 3 ; l++ )
					{
						double val = sx[ l ] * sy[ m ];
						xx += val * p[ l ].x;
						yy += val * p[ l ].y;
					}
				}

				xx *= _1_ax;
				yy *= _1_ay;

				double ct = -2000.0;
				if( xx < 0 || source.width( ) <= xx + 1 || yy < 0 || source.height( ) <= yy + 1 )
				{
				}
				else
				{
					difference_type ixx = static_cast< size_type >( xx );
					difference_type iyy = static_cast< size_type >( yy );
					typename source_image_type::const_pointer p = &source( ixx, iyy, 0 );
					xx -= ixx;
					yy -= iyy;

					ct = ( p[ d0 ] + ( p[ d3 ] - p[ d0 ] ) * xx ) + ( p[ d1 ] - p[ d0 ] + ( p[ d0 ] - p[ d1 ] + p[ d2 ] - p[ d3 ] ) * xx ) * yy;
				}
				target( x, y ) = static_cast< typename target_image_type::value_type >( ct );
			}
		}
	}


	// _pāC\[X摜^[Qbg摜ɕό`
	template < class TARGETTYPE, class TARGETTYPEA, class SOURCETYPE, class SOURCETYPEA, class CONTROLMESH >
	static void non_rigid_transformation_( array3< TARGETTYPE, TARGETTYPEA > &target, const array3< SOURCETYPE, SOURCETYPEA > &source, const CONTROLMESH &control_mesh,
											typename CONTROLMESH::difference_type mx = -1,
											typename CONTROLMESH::difference_type my = -1,
											typename CONTROLMESH::difference_type mz = -1,
											typename CONTROLMESH::size_type thread_id = 0,
											typename CONTROLMESH::size_type thread_num = 1
										)
	{
		typedef array3< TARGETTYPE, TARGETTYPEA >  target_image_type;
		typedef array3< SOURCETYPE, SOURCETYPEA >  source_image_type;
		typedef CONTROLMESH control_mesh_type;

		typedef typename target_image_type::size_type size_type;
		typedef typename target_image_type::difference_type difference_type;

		difference_type i, j, k;
		double u, v, w;
		double sx[ 4 ], sy[ 4 ], sz[ 4 ];
		double stepW = control_mesh.width( )  == 1 ? 1.0 : target.width( ) / static_cast< double >( control_mesh.width( ) - 1 );
		double stepH = control_mesh.height( ) == 1 ? 1.0 : target.height( ) / static_cast< double >( control_mesh.height( ) - 1 );
		double stepD = control_mesh.depth( )  == 1 ? 1.0 : target.depth( ) / static_cast< double >( control_mesh.depth( ) - 1 );
		double _1_stepW = 1.0 / stepW;
		double _1_stepH = 1.0 / stepH;
		double _1_stepD = 1.0 / stepH;
		double _1_ax = 1.0 / source.reso1( );
		double _1_ay = 1.0 / source.reso2( );
		double _1_az = 1.0 / source.reso2( );

		//stepW *= target.reso1( );
		//stepH *= target.reso2( );<commit_msg>Bスプライン変換に混入したバグを修正<commit_after>/// @file mist/registration.h
//!
//! @brief ́E񍄑̃WXg[V
//!
#ifndef __INCLUDE_MIST_REGISTRATION__
#define __INCLUDE_MIST_REGISTRATION__


#ifndef __INCLUDE_MIST_CONF_H__
#include "config/mist_conf.h"
#endif


#ifndef __INCLUDE_MIST_H__
#include "mist.h"
#endif

#ifndef __INCLUDE_MIST_MINIMIZATION__
#include "minimization.h"
#endif

#ifndef __INCLUDE_MIST_DRAWING__
#include "drawing.h"
#endif


// mistOԂ̎n܂
_MIST_BEGIN

namespace __non_rigid_registration_utility__
{
	inline void FFD( double v, double &f0, double &f1, double &f2, double &f3 )
	{
		double v_ = 1.0 - v;
		double v2 = v * v;
		double v3 = v2 * v;
		f0 = v_ * v_ * v_ / 6.0;
		f1 = 0.5 * v3 - v2 + 4.0 / 6.0;
		f2 = ( v2 - v3 + v ) * 0.5 + 1.0 / 6.0;
		f3 = v3 / 6.0;
	}


	// _pāC\[X摜^[Qbg摜ɕό`
	template < class TARGETTYPE, class TARGETTYPEA, class SOURCETYPE, class SOURCETYPEA, class CONTROLMESH >
	static void non_rigid_transformation_( array2< TARGETTYPE, TARGETTYPEA > &target, const array2< SOURCETYPE, SOURCETYPEA > &source, const CONTROLMESH &control_mesh,
											typename CONTROLMESH::difference_type mx = -1,
											typename CONTROLMESH::difference_type my = -1,
											typename CONTROLMESH::difference_type mz = -1,
											typename CONTROLMESH::size_type thread_id = 0,
											typename CONTROLMESH::size_type thread_num = 1
										)
	{
		typedef array2< TARGETTYPE, TARGETTYPEA >  target_image_type;
		typedef array2< SOURCETYPE, SOURCETYPEA >  source_image_type;
		typedef CONTROLMESH control_mesh_type;

		typedef typename target_image_type::size_type size_type;
		typedef typename target_image_type::difference_type difference_type;

		difference_type i, j;
		double u, v;
		double sx[ 4 ], sy[ 4 ];
		double stepW = control_mesh.width( ) == 1 ? 1.0 : target.width( ) / static_cast< double >( control_mesh.width( ) - 1 );
		double stepH = control_mesh.height( ) == 1 ? 1.0 : target.height( ) / static_cast< double >( control_mesh.height( ) - 1 );
		double _1_stepW = 1.0 / stepW;
		double _1_stepH = 1.0 / stepH;
		double _1_ax = 1.0 / source.reso1( );
		double _1_ay = 1.0 / source.reso2( );

		//stepW *= target.reso1( );
		//stepH *= target.reso2( );

		difference_type d0, d1, d2, d3;
		{
			difference_type cx = source.width( ) / 2;
			difference_type cy = source.height( ) / 2;
			typename source_image_type::const_pointer ppp = &source( cx, cy );
			d0 = 0;
			d1 = &source( cx    , cy + 1 ) - ppp;
			d2 = &source( cx + 1, cy + 1 ) - ppp;
			d3 = &source( cx + 1, cy     ) - ppp;
		}

		difference_type isx, isy, iex, iey;
		difference_type mw = control_mesh.width( );
		difference_type mh = control_mesh.height( );
		if( mx < 0 || mw <= mx || my < 0 || mh <= my )
		{
			isx = 0;
			iex = target.width( ) - 1;
			isy = 0;
			iey = target.height( ) - 1;
		}
		else
		{
			isx = static_cast< difference_type >( ( mx - 2 ) * stepW );
			iex = static_cast< difference_type >( ( mx + 2 ) * stepW );
			isy = static_cast< difference_type >( ( my - 2 ) * stepH );
			iey = static_cast< difference_type >( ( my + 2 ) * stepH );

			difference_type tw = target.width( );
			difference_type th = target.height( );
			isx = isx > 0  ? isx : 0;
			isy = isy > 0  ? isy : 0;
			iex = iex < tw ? iex : tw - 1;
			iey = iey < th ? iey : th - 1;
		}


		for( difference_type y = isy + thread_id ; y <= iey ; y += thread_num )
		{
			v = y * _1_stepH;
			j = static_cast< difference_type >( v );
			v -= j;
			j--;

			FFD( v, sy[ 0 ], sy[ 1 ], sy[ 2 ], sy[ 3 ] );

			for( difference_type x = isx ; x <= iex ; x++ )
			{
				u = x * _1_stepW;
				i = static_cast< difference_type >( u );
				u -= i;
				i--;

				FFD( u, sx[ 0 ], sx[ 1 ], sx[ 2 ], sx[ 3 ] );

				double xx = 0.0, yy = 0.0;
				for( difference_type m = 0 ; m <= 3 ; m++ )
				{
					typename control_mesh_type::const_pointer p = &control_mesh( i, j + m, 0 );
					for( difference_type l = 0 ; l <= 3 ; l++ )
					{
						double val = sx[ l ] * sy[ m ];
						xx += val * p[ l ].x;
						yy += val * p[ l ].y;
					}
				}

				xx *= _1_ax;
				yy *= _1_ay;

				double ct = -2000.0;
				if( xx < 0 || source.width( ) <= xx + 1 || yy < 0 || source.height( ) <= yy + 1 )
				{
				}
				else
				{
					difference_type ixx = static_cast< size_type >( xx );
					difference_type iyy = static_cast< size_type >( yy );
					typename source_image_type::const_pointer p = &source( ixx, iyy, 0 );
					xx -= ixx;
					yy -= iyy;

					ct = ( p[ d0 ] + ( p[ d3 ] - p[ d0 ] ) * xx ) + ( p[ d1 ] - p[ d0 ] + ( p[ d0 ] - p[ d1 ] + p[ d2 ] - p[ d3 ] ) * xx ) * yy;
				}
				target( x, y ) = static_cast< typename target_image_type::value_type >( ct );
			}
		}
	}


	// _pāC\[X摜^[Qbg摜ɕό`
	template < class TARGETTYPE, class TARGETTYPEA, class SOURCETYPE, class SOURCETYPEA, class CONTROLMESH >
	static void non_rigid_transformation_( array3< TARGETTYPE, TARGETTYPEA > &target, const array3< SOURCETYPE, SOURCETYPEA > &source, const CONTROLMESH &control_mesh,
											typename CONTROLMESH::difference_type mx = -1,
											typename CONTROLMESH::difference_type my = -1,
											typename CONTROLMESH::difference_type mz = -1,
											typename CONTROLMESH::size_type thread_id = 0,
											typename CONTROLMESH::size_type thread_num = 1
										)
	{
		typedef array3< TARGETTYPE, TARGETTYPEA >  target_image_type;
		typedef array3< SOURCETYPE, SOURCETYPEA >  source_image_type;
		typedef CONTROLMESH control_mesh_type;

		typedef typename target_image_type::size_type size_type;
		typedef typename target_image_type::difference_type difference_type;

		difference_type i, j, k;
		double u, v, w;
		double sx[ 4 ], sy[ 4 ], sz[ 4 ];
		double stepW = control_mesh.width( )  == 1 ? 1.0 : target.width( ) / static_cast< double >( control_mesh.width( ) - 1 );
		double stepH = control_mesh.height( ) == 1 ? 1.0 : target.height( ) / static_cast< double >( control_mesh.height( ) - 1 );
		double stepD = control_mesh.depth( )  == 1 ? 1.0 : target.depth( ) / static_cast< double >( control_mesh.depth( ) - 1 );
		double _1_stepW = 1.0 / stepW;
		double _1_stepH = 1.0 / stepH;
		double _1_stepD = 1.0 / stepD;
		double _1_ax = 1.0 / source.reso1( );
		double _1_ay = 1.0 / source.reso2( );
		double _1_az = 1.0 / source.reso3( );

		//stepW *= target.reso1( );
		//stepH *= target.reso2( );
{'fraction_non_alphanumeric': 0.1415717811927701, 'fraction_numerical': 0.02977944103829923, 'mean_word_length': 2.7039221468593335, 'pattern_counts': {'":': 0, '<': 79, '<?xml version=': 0, '>': 55, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3175397', 'n_tokens_mistral': 6056, 'n_tokens_neox': 5006, 'n_words': 1260}