/*---------------------------------------------------------------------------
   
   File name 	:	oddeven.c
   Author 		:	N. Devillard
   Created on	:	Oct 2001
   Description	:	Filter to remove odd-even effect on ISAAC

 *--------------------------------------------------------------------------*/

/*
	$Id: oe_filt.c,v 1.2 2002/02/08 10:09:18 ndevilla Exp $
	$Author: ndevilla $
	$Date: 2002/02/08 10:09:18 $
	$Revision: 1.2 $
*/

/*---------------------------------------------------------------------------
   								Includes
 ---------------------------------------------------------------------------*/

#include "image_handling.h"
#include "cube_handling.h"
#include "extraction.h"
#include "image_intops.h"
#include "fourier.h"



/*---------------------------------------------------------------------------
  							Function codes
 ---------------------------------------------------------------------------*/

/*-------------------------------------------------------------------------*/
/**
  @brief    Remove the odd-even effect inside an image quadrant.
  @param    im  Image to process
  @return   1 newly allocated image.

  This function tries to remove the odd-even effect from a single
  quadrant. The input image size must be a power of 2 and the image
  must be square.

  THIS IS HIGHLY ISAAC-SPECIFIC!!
 */
/*--------------------------------------------------------------------------*/
static image_t * image_de_oddeven_quad(image_t * im)
{
	image_t	*	cleaned ;
	cube_t	*	freq_i ;
	cube_t	*	freq_o ;
	cube_t	*	freq_i_amp ;
	int			lx, ly ;

	if (im==NULL) return NULL ;

	lx = im->lx ;
	ly = im->ly ;

	/* Apply FFT to input image */
	freq_i = image_fft(im, NULL, FFT_FORWARD);
	/* Convert to amplitude/phase */
	freq_i_amp = cube_conv_xy_rtheta(freq_i);
	cube_del(freq_i);
	/* Nullify the odd-even frequencies */
	freq_i_amp->plane[0]->data[lx/2] = 0 ;
	freq_i_amp->plane[0]->data[lx/2 + (ly-1)*lx] = 0 ;
	/* Convert to X/Y */
	freq_i = cube_conv_rtheta_xy(freq_i_amp);
	cube_del(freq_i_amp);
	/* FFT back to image space */
	freq_o = image_fft(freq_i->plane[0], freq_i->plane[1], FFT_INVERSE);
	cube_del(freq_i);
	cleaned = freq_o->plane[0] ;
	freq_o->plane[0] = NULL ;
	cube_del(freq_o);
	return cleaned ;
}

/*-------------------------------------------------------------------------*/
/**
  @brief    Remove the odd-even effect inside an ISAAC image.
  @param    im  Image to process
  @return   1 newly allocated image.

  This function tries to remove the odd-even effect from an ISAAC
  image, applying image_de_oddeven_quad() to each quadrant in
  sequence. The input image size must be a power of 2 and the image
  must be square.

  THIS IS HIGHLY ISAAC-SPECIFIC!!
 */
/*--------------------------------------------------------------------------*/
image_t * image_de_oddeven(image_t * im)
{
	image_t	*	quad ;
	image_t	*	f_quad ;
	image_t	*	cleaned ;
	image_t	*	cleaned_tmp ;

	cleaned = image_new(im->lx, im->ly);

	compute_status("filtering odd-even effect", 0, 4, 1);
	quad = image_getvig(im, 1, 1, im->lx/2, im->ly/2);
	f_quad = image_de_oddeven_quad(quad);
	image_del(quad);
	cleaned_tmp =
	image_paste_vig(cleaned, f_quad, 1, 1,
					1, 1, f_quad->lx, f_quad->ly);
	image_del(f_quad);
	image_del(cleaned) ;
	cleaned = cleaned_tmp ;

	compute_status("filtering odd-even effect", 1, 4, 1);
	quad = image_getvig(im, 1+im->lx/2, 1, im->lx, im->ly/2);
	f_quad = image_de_oddeven_quad(quad);
	image_del(quad);
	cleaned_tmp =
	image_paste_vig(cleaned, f_quad, 1+im->lx/2, 1,
					1, 1, f_quad->lx, f_quad->ly);
	image_del(f_quad);
	image_del(cleaned);
	cleaned = cleaned_tmp ;

	compute_status("filtering odd-even effect", 2, 4, 1);
	quad = image_getvig(im, 1, 1+im->ly/2, im->lx/2, im->ly);
	f_quad = image_de_oddeven_quad(quad);
	image_del(quad);
	cleaned_tmp =
	image_paste_vig(cleaned, f_quad, 1, 1+im->ly/2,
					1, 1, f_quad->lx, f_quad->ly);
	image_del(f_quad);
	image_del(cleaned);
	cleaned = cleaned_tmp ;

	compute_status("filtering odd-even effect", 3, 4, 1);
	quad = image_getvig(im, 1+im->lx/2, 1+im->ly/2, im->lx, im->ly);
	f_quad = image_de_oddeven_quad(quad);
	image_del(quad);
	cleaned_tmp =
	image_paste_vig(cleaned, f_quad, 1+im->lx/2, 1+im->ly/2,
					1, 1, f_quad->lx, f_quad->ly);
	image_del(f_quad);
	image_del(cleaned);
	cleaned = cleaned_tmp ;

	return cleaned ;

}