stereoview::image_id_type stereoview::add_image(const image_type& img, const option<float>& focal_length)
{
    if (!stored_image_db.empty() && (img.width() != width() || img.height() != height()))
        throw localized_invalid_argument(HERE(nfmt<4>("image size is %1 x %2 while expected size is %3 x %4") (img.width()) (img.height()) (width()) (height())));
    const unsigned int id = fresh_int();
    add_stored_image(id, img, focal_length);
    return id;
}

IonDetector::result_type IonDetector::blob_detector( const image_type& img ) {
	Kernel logs[3];
	float sigmas[3];

	logs[0] = Kernel::LoG( 4, 0.6 );
	sigmas[0] = 0.85;

	logs[1] = Kernel::LoG( 6, 1.0 );
	sigmas[1] = 1.0;

	logs[2] = Kernel::LoG( 20, 5.0 );
	sigmas[2] = 2.;

	image_type img_result( img.extents );
	result_type results;

	const size_t hsize = 2000;
	std::vector< size_t > histogram( hsize, 0 );

	size_t g = 2;
	logs[g].apply_kernel( img, img_result );
	save_file( "log.fits", img_result );

	float img_max = *std::max_element( img_result.ptr(), 
		img_result.ptr() + img.num_elements() );
	for( 	float* iter = img_result.ptr(); 
			iter != img_result.ptr() + img.num_elements(); 
			++iter ) {

		int p = (size_t)( ( *iter + img_max ) / 2 / img_max * hsize ); 
		if( p >= (int)hsize ) 	p = hsize-1;
		if( p < 0 )				p = 0;
		histogram[ p ] ++;
	}

	size_t nelems = 0, counter = hsize-1;
	while( nelems < (size_t)blob_threshold )
		nelems += histogram[ counter-- ];

	float thresh = (float)counter / hsize * img_max - img_max / 2.;
	for( image_type::index i = 12; i < img.extents.first-12; ++i ) 
	for( image_type::index j = 12; j < img.extents.second-12; ++j ) {
		if( img_result( i, j ) > thresh &&
			img_result( i, j ) > img_result( i+1, j ) &&
			img_result( i, j ) > img_result( i-1, j ) &&
			img_result( i, j ) > img_result( i, j+1 ) &&
			img_result( i, j ) > img_result( i, j-1 ) &&
			img_result( i, j ) > img_result( i+1, j+1 ) &&
			img_result( i, j ) > img_result( i+1, j-1 ) &&
			img_result( i, j ) > img_result( i-1, j+1 ) &&
			img_result( i, j ) > img_result( i-1, j-1 ) ) {

			results.push_back( IonData( i, j, sigmas[g] ) );
		}
	}

	return results;
}

 unsigned compare(image_type const & actual, boost::filesystem::path const& reference) const
 {
     std::ifstream stream(reference.string().c_str(), std::ios_base::in | std::ios_base::binary);
     if (!stream)
     {
         throw std::runtime_error("Could not open: " + reference.string());
     }
     std::string expected(std::istreambuf_iterator<char>(stream.rdbuf()), std::istreambuf_iterator<char>());
     return std::max(actual.size(), expected.size()) - std::min(actual.size(), expected.size());
 }

void fft_round_up(image_type& I,pos_type& from,pos_type& to)
{
    image_type newI(fft_round_up_geometry(I.geometry()));
    for(int dim = 0;dim < image_type::dimension;++dim)
    {
        from[dim] = (newI.geometry()[dim]-I.geometry()[dim]) >> 1;
        to[dim] = from[dim] + I.geometry()[dim];
    }
    image::draw(I,newI,from);
    I.swap(newI);
}

void Image2DIB(const image_type &image, BYTE *dib)
{
	using namespace rss;
	const int all_header_size = 14 + 40 +1024;
	size_t dib_size = all_header_size 
		+ ( image.width() + ( (image.width()%4)?(4-image.width()%4):0) )* image.height();	
	
	std::ostrstream output(reinterpret_cast<char *>(dib), dib_size);
	BMPImageIO<image_type> image_io;
	if(!image_io.write(output, image))
		throw rss::Exception("library can not write this dib");
}

IonDetector::result_type IonDetector::hough_transform( image_type& edges ) {
	const unsigned short num_radii = 4;

	std::vector< unsigned short > hough( 
		edges.num_elements()*num_radii, 0 );

	for( size_t r = 2; r < 6; ++r ) 
	for( image_type::index i = r+1; i<edges.extents.first-r-1; ++i ) 
	for( image_type::index j = r+1; j<edges.extents.second-r-1; ++j ) {
		if( edges( i, j ) ) {
			for( float theta = 0.; theta < 6.3; theta += 1. / r ) {
				size_t x = (size_t)( i + (float)r*cos( theta ) + 0.5f);
				size_t y = (size_t)( j + (float)r*sin( theta ) + 0.5f);
				hough[ (r-2)*edges.num_elements() +
					y*edges.extents.first + x ]++;
			}
		}
	}
		
	result_type results;
	for( size_t r = 5; r >= 2; --r )
	for( image_type::index i = 8; i < edges.extents.first-8; ++i ) 
	for( image_type::index j = 8; j < edges.extents.second-8; ++j ) {
		size_t index = (r-2)*edges.num_elements() + 
			j*edges.extents.first + i;
		const unsigned short value = hough[index];
		if( value > (size_t)(hough_threshold * r) ) {
			size_t h = 0;
			for( ;  h < results.size(); ++h ) {
				size_t xsq = (i - results[h].x)*(i - results[h].x);
				size_t ysq = (j - results[h].y)*(j - results[h].y);
				if( xsq + ysq < 36 ) break;
			}
			
			if( h != results.size() ) continue;
			if( hough[index + edges.extents.first] > value ) continue;
			if( hough[index - edges.extents.first] > value ) continue;
			if( hough[index - 1] > value ) continue;
			if( hough[index + 1] > value ) continue;

			results.push_back( IonData( i, j, r ) );
			std::cout<< "Circle " << i << " " << j << " " << r << std::endl;
			std::cout<< value << std::endl;
		}
	}
	
	return results;
}

    void convert(mapnik::grid::data_type const & grid, image_type & image) const
    {
        for (std::size_t y = 0; y < grid.height(); ++y)
        {
            mapnik::grid::value_type const * grid_row = grid.get_row(y);
            image_type::pixel_type * image_row = image.get_row(y);
            for (std::size_t x = 0; x < grid.width(); ++x)
            {
                mapnik::grid::value_type val = grid_row[x];

                if (val == mapnik::grid::base_mask)
                {
                    image_row[x] = 0;
                    continue;
                }
                if (val < 0)
                {
                    throw std::runtime_error("grid renderer: feature id is negative.");
                }

                val *= 100;

                if (val > 0x00ffffff)
                {
                    throw std::runtime_error("grid renderer: feature id is too high.");
                }

                image_row[x] = val | 0xff000000;
            }
        }
    }

void RSS_MultiSensor_PCA_Fusion(const image_type &image1, const image_type &image2, image_type &result)
{
	using namespace rss;
	rss::ImageVector<image_type> input_vector;
	input_vector.push_back(image1);
	if(image1.size() != image2.size()) {
		BilinearInterpolation<image_type> interpolate(image1.size());
		image_type temp;
		interpolate(image2, temp);
		input_vector.push_back(temp);
	} else { 
		input_vector.push_back(image2);
	}

	PCAFusion<image_type> fusion(image1.width(), image1.height());
	fusion(input_vector, result);		
}

void RSS_MultiSensor_Wavelet_Fusion(const image_type &image1, const image_type &image2, image_type &result)
{
	using namespace rss;
	rss::ImageVector<image_type> input_vector;
	input_vector.push_back(image1);
	if(image1.size() != image2.size()) {
		BilinearInterpolation<image_type> interpolate(image1.size());
		image_type result;
		interpolate(image2, result);
		input_vector.push_back(result);
	} else { 
		input_vector.push_back(image2);
	}

	WaveletFusion<image_type> fusion(image1.size(), FilterSet::Haar());
	fusion(input_vector, result);		

}

void RSS_MultiSensor_Contour_Register(const image_type &VLImage, const image_type &IRImage, float para[6], image_type &result)
{
	using namespace rss;
	ContourRegister contour_register(HomoModel::Similarity, 3.1, 0.0, 0.0, 0.0, 30, 0.18, 0.8, 0.5);
	HomoModel model; 
	contour_register(VLImage, IRImage,model);
	HomoTrans<image_type> homo_trans(model, VLImage.size());
	homo_trans.inverse(IRImage, result);
}

void fft_shift_x(image_type& I)
{
    int half_w = I.width() >> 1;
    int half_w_1 = half_w-1;
    int w_1 = I.width()-1;
    int quater_w = half_w >> 1;
    typename image_type::iterator iter1 = I.begin();
    typename image_type::iterator end = I.end();
    for(;iter1 != end;iter1 += I.width())
    {
        typename image_type::iterator iter2 = iter1+half_w;
        for(int x = 0,rx = half_w_1;x < quater_w;++x,--rx)
        {
            std::swap(iter1[x],iter1[rx]);
            std::swap(iter2[x],iter2[rx]);
        }
    }
}

void fft_shift_z(image_type& I)
{
    int wh = I.plane_size();
    int half_size = I.size() >> 1;
    int half_size_1 = half_size-wh;
    int quater_size = half_size >> 1;
    typename image_type::iterator iter1 = I.begin();
    typename image_type::iterator iter2 = iter1+half_size;
    typename image_type::iterator end = iter1+wh;
    for(;iter1 != end;++iter1,++iter2)
    {
        for(int z = 0,rz = half_size_1;z < quater_size;z+=wh,rz-=wh)
        {
            std::swap(iter1[z],iter1[rz]);
            std::swap(iter2[z],iter2[rz]);
        }
    }
}

 unsigned compare(image_type const & actual, boost::filesystem::path const& reference) const
 {
     std::ifstream stream(reference.string().c_str(),std::ios_base::in|std::ios_base::binary);
     if (!stream.is_open())
     {
         throw std::runtime_error("could not open: '" + reference.string() + "'");
     }
     std::string expected(std::istreambuf_iterator<char>(stream.rdbuf()),(std::istreambuf_iterator<char>()));
     stream.close();
     return std::fabs(actual.size() - expected.size());
 }

void decode_image_impl(image_type& image) {
  if (image.m_format == Format::RAW_ARRAY) {
    return;
  }

  char* buf = NULL;
  size_t length = 0;

  if (image.m_format == Format::JPG) {
    decode_jpeg((const char*)image.get_image_data(), image.m_image_data_size,
                &buf, length);
  } else if (image.m_format == Format::PNG) {
    decode_png((const char*)image.get_image_data(), image.m_image_data_size,
                &buf, length);
  } else {
    log_and_throw(std::string("Cannot decode image. Unknown format."));
  };
  image.m_image_data.reset(buf);
  image.m_image_data_size = length;
  image.m_format = Format::RAW_ARRAY;
}

void save_file( const std::string& filename, const image_type& img ) {
	fitsfile* f;
	int status = 0;

    std::string fname = "!" + filename;
	fits_create_file( &f, fname.c_str(), &status );
	
	long dims[2];
	dims[0] = img.extents.first;	dims[1] = img.extents.second;

	typedef CFitsIOTypeMap<	image_type::element > map;
	fits_create_img( f, map::image_type, 2, dims, &status );
	fits_write_img( f, map::type, 1, img.num_elements(), 
		(void*)img.ptr(), &status );
        
    fits_close_file( f, &status );

	if( status ) {
		fits_report_error( stderr, status );
	}
}

void fft_shift_y(image_type& I)
{
    int w = I.width();
    int half_wh = I.plane_size() >> 1;
    int half_wh_1 = half_wh-w;
    int quater_wh = half_wh >> 1;
    typename image_type::iterator iter1 = I.begin();
    typename image_type::iterator end = I.end();
    for(;iter1 != end;iter1 += I.plane_size())
    {
        typename image_type::iterator iter1_x = iter1;
        typename image_type::iterator iter1_x_end = iter1+w;
        typename image_type::iterator iter2_x = iter1_x+half_wh;
        for(;iter1_x != iter1_x_end;++iter1_x,++iter2_x)
        for(int y = 0,ry = half_wh_1;y < quater_wh;y+=w,ry-=w)
        {
            std::swap(iter1_x[y],iter1_x[ry]);
            std::swap(iter2_x[y],iter2_x[ry]);
        }
    }
}

void encode_image_impl(image_type& image) {
  if (image.m_format != Format::RAW_ARRAY){
    return;
  } 

  char* buf = NULL;
  size_t length = 0;

  encode_png((const char*)image.get_image_data(), image.m_width, image.m_height, image.m_channels, &buf, length);
  image.m_image_data.reset(buf);
  image.m_image_data_size = length;
  image.m_format = Format::PNG;
}

    inline void load (
        const image_type& img
    )
    {
        feat_image.set_size(img.nr(), img.nc());
        assign_all_pixels(feat_image,0);
        for (long r = 1; r+1 < img.nr(); ++r)
        {
            for (long c = 1; c+1 < img.nc(); ++c)
            {
                unsigned char f = 0;
                if (img[r][c])   f |= 0x1;
                if (img[r][c+1]) f |= 0x2;
                if (img[r][c-1]) f |= 0x4;
                if (img[r+1][c]) f |= 0x8;
                if (img[r-1][c]) f |= 0x10;

                // Store the code value for the pattern of pixel values in the 4-connected
                // neighborhood around this row and column.
                feat_image[r][c] = f;
            }
        }
    }

    void operator()(image_type& src)
    {
        typedef tipl::image<typename image_type::value_type,image_type::dimension> image_buf_type;
        image_buf_type gx;
        gradient_2x(src,gx);

        image_buf_type gy;
        gradient_2y(src,gy);

        for(size_t index = 0;index < src.size();++index)
        {
            float fx = gx[index];
            float fy = gy[index];
            src[index] = std::sqrt(fx*fx+fy*fy);
        }
    }

void RSS_MultiSensor_FFT_Register(const image_type &VLImage, const image_type &IRImage, float para[6], image_type &result)
{
	using namespace rss;
	image_type temp;
	if(para) {
		double scale1 = para[0] / para[3] * para[4] / para[1] * VLImage.width() / IRImage.width();
		double scale2 = para[0] / para[3] * para[5] / para[2] * VLImage.height() / IRImage.height();
		double scale = sqrt(scale1*scale2);

		HomoModel pre_trans_model;
		pre_trans_model.SetSimilarity(scale, 0, IRImage.width() * (1 - scale) / 2.0, IRImage.height() * (1 - scale) / 2.0);
		HomoTrans<image_type> pre_trans(pre_trans_model, VLImage.size());
		pre_trans.operator ()(IRImage, temp);
	} else {
		temp = IRImage;
	}
	HomoModel model;
	SimilarityEstimation(true, true, false, SimilarityEstimation::OP_NONE, SimilarityEstimation::FILTER_NONE)(VLImage, temp, model);
	HomoTrans<image_type> homo_trans(model, VLImage.size());
	homo_trans.operator()(temp, result);

}