/*
 * itrans-zaefferer.c
 *
 * Zaefferer peak search
 *
 * (c) 2007 Thomas White <taw27@cam.ac.uk>
 *
 *  dtr - Diffraction Tomography Reconstruction
 *
 */

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

#include <stdint.h>

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

#define PEAK_WINDOW_SIZE 20

unsigned int itrans_peaksearch_zaefferer(int16_t *image, ControlContext *ctx, double tilt_degrees, ImageDisplay *imagedisplay) {

	int x, y;
	unsigned int n_reflections;
	
	n_reflections = 0;
	for ( x=1; x<ctx->width-1; x++ ) {
		for ( y=1; y<ctx->height-1; y++ ) {

			double dx1, dx2, dy1, dy2;
			double dxs, dys;
			double grad;

			if ( ctx->max_d ) {
				if ( (x-ctx->x_centre)*(x-ctx->x_centre) + (y-ctx->y_centre)*(y-ctx->y_centre) > ctx->max_d*ctx->max_d ) {
					continue;
				}
			}
			
			/* Get gradients */
			dx1 = image[x+ctx->width*y] - image[(x+1)+ctx->width*y];
			dx2 = image[(x-1)+ctx->width*y] - image[x+ctx->width*y];
			dy1 = image[x+ctx->width*y] - image[(x+1)+ctx->width*(y+1)];
			dy2 = image[x+ctx->width*(y-1)] - image[x+ctx->width*y];
			
			/* Average gradient measurements from both sides */
			dxs = ((dx1*dx1) + (dx2*dx2)) / 2;
			dys = ((dy1*dy1) + (dy2*dy2)) / 2;
			
			/* Calculate overall gradient */
			grad = dxs + dys;

			if ( grad > 400 ) {
			
				unsigned int mask_x, mask_y;
				unsigned int sx, sy;
				double max;
				unsigned int did_something = 1;
				
				mask_x = x;
				mask_y = y;
				
				while ( (did_something) && (distance(mask_x, mask_y, x, y)<50) ) {
					max = image[mask_x+ctx->width*mask_y];
					did_something = 0;
					for ( sy=biggest(mask_y-PEAK_WINDOW_SIZE/2, 0); sy<smallest(mask_y+PEAK_WINDOW_SIZE/2, ctx->height); sy++ ) {
						for ( sx=biggest(mask_x-PEAK_WINDOW_SIZE/2, 0); sx<smallest(mask_x+PEAK_WINDOW_SIZE/2, ctx->width); sx++ ) {
							if ( image[sx+ctx->width*sy] > max ) {
								max = image[sx+ctx->width*sy];
								mask_x = sx;
								mask_y = sy;
								did_something = 1;
							}
						}
					}
				}
				
				if ( !did_something ) {
					if ( ctx->fmode == FORMULATION_PIXELSIZE ) {
						reflection_add_from_reciprocal(ctx, (signed)(mask_x-ctx->width/2)*ctx->pixel_size,
										(signed)(mask_y-ctx->height/2)*ctx->pixel_size,
										tilt_degrees, image[mask_x + ctx->width*mask_y]);
					} else {
						reflection_add_from_dp(ctx, (signed)(mask_x-ctx->width/2), (signed)(mask_y-ctx->height/2),
										tilt_degrees, image[mask_x + ctx->width*mask_y]);
					}
					if ( ctx->first_image ) {
						imagedisplay_mark_point(imagedisplay, mask_x, mask_y);
					}
					n_reflections++;
				}
				
			}
		}
	}

	return n_reflections;

}