/*
 * itrans-threshold.c
 *
 * Threshold and adaptive threshold peak searches
 *
 * (c) 2007 Thomas White <taw27@cam.ac.uk>
 *
 *  dtr - Diffraction Tomography Reconstruction
 *
 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdint.h>
#include <assert.h>

#include "control.h"
#include "imagedisplay.h"
#include "reflections.h"

unsigned int itrans_peaksearch_threshold(ImageRecord imagerecord, ControlContext *ctx) {

	int width, height;
	int x, y;
	unsigned int n_reflections = 0;
	uint16_t *image = imagerecord.image;
	double tilt_degrees = imagerecord.tilt;
	uint16_t max = 0;
	
	width = imagerecord.width;
	height = imagerecord.height;
	
	/* Simple Thresholding */
	for ( y=0; y<height; y++ ) {
		for ( x=0; x<width; x++ ) {
			if ( image[x + width*y] > max ) max = image[x + width*y];
		}
	}
	
	for ( y=0; y<height; y++ ) {
		for ( x=0; x<width; x++ ) {
			if ( image[x + width*y] > max/3 ) {
				assert(x<width);
				assert(y<height);
				assert(x>=0);
				assert(y>=0);
				if ( imagerecord.fmode == FORMULATION_PIXELSIZE ) {
					reflection_add_from_reciprocal(ctx,
								(signed)(x-imagerecord.x_centre)*imagerecord.pixel_size,
								(signed)(y-imagerecord.y_centre)*imagerecord.pixel_size,
								tilt_degrees, image[x + width*y]);
				} else {
					reflection_add_from_dp(ctx,
								(signed)(x-imagerecord.x_centre)/imagerecord.resolution,
								(signed)(y-imagerecord.y_centre)/imagerecord.resolution,
								tilt_degrees, image[x + width*y]);
				}
				n_reflections++;
			}
		}
	}
	
	return n_reflections;
	
}

unsigned int itrans_peaksearch_adaptive_threshold(ImageRecord imagerecord, ControlContext *ctx) {

	uint16_t max_val = 0;
	int width, height;
	unsigned int n_reflections = 0;
	uint16_t *image;
	double tilt_degrees = imagerecord.tilt;
	uint16_t max;
	int x, y;
	
	image = imagerecord.image;
	width = imagerecord.width;
	height = imagerecord.height;
	
	max = 0;
	for ( y=0; y<height; y++ ) {
		for ( x=0; x<width; x++ ) {
			if ( image[x + width*y] > max ) max = image[x + width*y];
		}
	}
	
	/* Adaptive Thresholding */
	do {
	
		int max_x = 0;
		int max_y = 0;;
		
		/* Locate the highest point */
		max_val = 0;
		for ( y=0; y<height; y++ ) {
			for ( x=0; x<width; x++ ) {

				if ( image[x + width*y] > max_val ) {
					max_val = image[x + width*y];
					max_x = x;  max_y = y;
				}
			
			}
		}
		
		if ( max_val > max/10 ) {
			assert(max_x<width);
			assert(max_y<height);
			assert(max_x>=0);
			assert(max_y>=0);
			if ( imagerecord.fmode == FORMULATION_PIXELSIZE ) {
				reflection_add_from_reciprocal(ctx,
								(signed)(max_x-imagerecord.x_centre)*imagerecord.pixel_size,
								(signed)(max_y-imagerecord.y_centre)*imagerecord.pixel_size,
								tilt_degrees, image[max_x + width*max_y]);
			} else {
				reflection_add_from_dp(ctx, (signed)(max_x-imagerecord.x_centre)/imagerecord.resolution,
							(signed)(max_y-imagerecord.y_centre)/imagerecord.resolution,
							tilt_degrees, image[max_x + width*max_y]);
			}
			n_reflections++;
			
			/* Remove it and its surroundings */
			for ( y=((max_y-10>0)?max_y-10:0); y<((max_y+10)<height?max_y+10:height); y++ ) {
				for ( x=((max_x-10>0)?max_x-10:0); x<((max_x+10)<width?max_x+10:width); x++ ) {
					image[x + width*y] = 0;
				}
			}
		}
		
	} while ( max_val > 50 );
	
	return n_reflections;

}