/* Inverse FFT

 *

 * Author: Nicos Dessipris

 * Written on: 12/04/1990

 * Modified on :

 * 28/6/95 JC

 *  - rewritten, based on new im_invfft() code

 * 10/9/98 JC

 *  - frees memory more quickly

 * 2/4/02 JC

 *  - fftw code added

 * 13/7/02 JC

 *  - Type reset

 * 27/2/03 JC

 *  - tiny speed-up ... save 1 copy on write

 * 22/1/04 JC

 *  - oops, fix for segv on wider than high fftw transforms

 * 3/11/04

 *  - added fftw3 support

 * 7/2/10

 *  - gtkdoc

 * 27/1/12

 *  - better setting of interpretation

 *  - remove own fft fallback code

 *  - remove fftw2 path

 *  - reduce memuse

 * 3/1/14

 *  - redone as a class

 * 15/12/23 [akash-akya]

 *  - add locks

 */

/*

  This file is part of VIPS.

  VIPS is free software; you can redistribute it and/or modify

  it under the terms of the GNU Lesser General Public License as published by

  the Free Software Foundation; either version 2 of the License, or

  (at your option) any later version.

  This program is distributed in the hope that it will be useful,

  but WITHOUT ANY WARRANTY; without even the implied warranty of

  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

  GNU Lesser General Public License for more details.

  You should have received a copy of the GNU Lesser General Public License

  along with this program; if not, write to the Free Software

  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA

  02110-1301  USA

 */

/*

  These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk

 */

#ifdef HAVE_CONFIG_H

#include <config.h>

#endif /*HAVE_CONFIG_H*/

#include <glib/gi18n-lib.h>

#include <stdio.h>

#include <stdlib.h>

#include <vips/vips.h>

#include <vips/internal.h>

#include "pfreqfilt.h"

#ifdef HAVE_FFTW

#include <fftw3.h>

typedef struct _VipsInvfft {

  VipsFreqfilt parent_instance;

  gboolean real;

} VipsInvfft;

typedef VipsFreqfiltClass VipsInvfftClass;

G_DEFINE_TYPE(VipsInvfft, vips_invfft, VIPS_TYPE_FREQFILT);

/* Complex to complex inverse transform.

 */

static int

cinvfft1(VipsObject *object, VipsImage *in, VipsImage **out)

{

  VipsImage **t = (VipsImage **) vips_object_local_array(object, 4);

  VipsInvfft *invfft = (VipsInvfft *) object;

  VipsObjectClass *class = VIPS_OBJECT_GET_CLASS(invfft);

  fftw_plan plan;

  double *planner_scratch;

  if (vips_check_mono(class->nickname, in) ||

    vips_check_uncoded(class->nickname, in))

    return -1;

  /* Convert input to a complex double membuffer.

   */

  *out = vips_image_new_memory();

  if (vips_cast_dpcomplex(in, &t[0], NULL) ||

    vips_image_write(t[0], *out))

    return -1;

  /* Make the plan for the transform. Yes, they really do use nx for

   * height and ny for width.

   */

  if (!(planner_scratch = VIPS_ARRAY(invfft,

        VIPS_IMAGE_N_PELS(in) * 2, double)))

    return -1;

  g_mutex_lock(&vips__fft_lock);

  if (!(plan = fftw_plan_dft_2d(in->Ysize, in->Xsize,

        (fftw_complex *) planner_scratch,

        (fftw_complex *) planner_scratch,

        FFTW_BACKWARD,

        0))) {

    g_mutex_unlock(&vips__fft_lock);

    vips_error(class->nickname,

      "%s", _("unable to create transform plan"));

    return -1;

  }

  g_mutex_unlock(&vips__fft_lock);

  fftw_execute_dft(plan,

    (fftw_complex *) (*out)->data, (fftw_complex *) (*out)->data);

  g_mutex_lock(&vips__fft_lock);

  fftw_destroy_plan(plan);

  g_mutex_unlock(&vips__fft_lock);

  (*out)->Type = VIPS_INTERPRETATION_B_W;

  return 0;

}

/* Complex to real inverse transform.

 */

static int

rinvfft1(VipsObject *object, VipsImage *in, VipsImage **out)

{

  VipsImage **t = (VipsImage **) vips_object_local_array(object, 4);

  VipsInvfft *invfft = (VipsInvfft *) object;

  VipsObjectClass *class = VIPS_OBJECT_GET_CLASS(invfft);

  const int half_width = in->Xsize / 2 + 1;

  double *half_complex;

  double *planner_scratch;

  fftw_plan plan;

  int x, y;

  double *q, *p;

  /* Convert input to a complex double membuffer.

   */

  t[1] = vips_image_new_memory();

  if (vips_cast_dpcomplex(in, &t[0], NULL) ||

    vips_image_write(t[0], t[1]))

    return -1;

  /* Build half-complex image.

   */

  if (!(half_complex = VIPS_ARRAY(invfft,

        t[1]->Ysize * half_width * 2, double)))

    return -1;

  q = half_complex;

  for (y = 0; y < t[1]->Ysize; y++) {

    p = ((double *) t[1]->data) + (guint64) y * t[1]->Xsize * 2;

    for (x = 0; x < half_width; x++) {

      q[0] = p[0];

      q[1] = p[1];

      p += 2;

      q += 2;

    }

  }

  /* Make mem buffer real image for output.

   */

  *out = vips_image_new_memory();

  if (vips_image_pipelinev(*out, VIPS_DEMAND_STYLE_ANY, t[1], NULL))

    return -1;

  (*out)->BandFmt = VIPS_FORMAT_DOUBLE;

  (*out)->Type = VIPS_INTERPRETATION_B_W;

  if (vips_image_write_prepare(*out))

    return -1;

  /* Make the plan for the transform. Yes, they really do use nx for

   * height and ny for width.

   */

  if (!(planner_scratch = VIPS_ARRAY(invfft,

        t[1]->Ysize * half_width * 2, double)))

    return -1;

  g_mutex_lock(&vips__fft_lock);

  if (!(plan = fftw_plan_dft_c2r_2d(t[1]->Ysize, t[1]->Xsize,

        (fftw_complex *) planner_scratch, (double *) (*out)->data,

        0))) {

    g_mutex_unlock(&vips__fft_lock);

    vips_error(class->nickname,

      "%s", _("unable to create transform plan"));

    return -1;

  }

  g_mutex_unlock(&vips__fft_lock);

  fftw_execute_dft_c2r(plan,

    (fftw_complex *) half_complex, (double *) (*out)->data);

  g_mutex_lock(&vips__fft_lock);

  fftw_destroy_plan(plan);

  g_mutex_unlock(&vips__fft_lock);

  return 0;

}

static int

vips_invfft_build(VipsObject *object)

{

  VipsFreqfilt *freqfilt = VIPS_FREQFILT(object);

  VipsInvfft *invfft = (VipsInvfft *) object;

  VipsImage **t = (VipsImage **) vips_object_local_array(object, 4);

  VipsImage *in;

  if (VIPS_OBJECT_CLASS(vips_invfft_parent_class)->build(object))

    return -1;

  in = freqfilt->in;

  if (vips_image_decode(in, &t[0]))

    return -1;

  in = t[0];

  if (invfft->real) {

    if (vips__fftproc(VIPS_OBJECT(invfft),

        in, &t[1], rinvfft1))

      return -1;

  }

  else {

    if (vips__fftproc(VIPS_OBJECT(invfft),

        in, &t[1], cinvfft1))

      return -1;

  }

  if (vips_image_write(t[1], freqfilt->out))

    return -1;

  return 0;

}

static void

vips_invfft_class_init(VipsInvfftClass *class)

{

  GObjectClass *gobject_class = G_OBJECT_CLASS(class);

  VipsObjectClass *vobject_class = VIPS_OBJECT_CLASS(class);

  gobject_class->set_property = vips_object_set_property;

  gobject_class->get_property = vips_object_get_property;

  vobject_class->nickname = "invfft";

  vobject_class->description = _("inverse FFT");

  vobject_class->build = vips_invfft_build;

  VIPS_ARG_BOOL(class, "real", 4,

    _("Real"),

    _("Output only the real part of the transform"),

    VIPS_ARGUMENT_OPTIONAL_INPUT,

    G_STRUCT_OFFSET(VipsInvfft, real),

    FALSE);

}

static void

vips_invfft_init(VipsInvfft *invfft)

{

}

#endif /*HAVE_FFTW*/

/**

 * vips_invfft: (method)

 * @in: input image

 * @out: (out): output image

 * @...: %NULL-terminated list of optional named arguments

 *

 * Optional arguments:

 *

 * * @real: only output the real part

 *

 * Transform an image from Fourier space to real space. The result is complex.

 * If you are OK with a real result, set @real, it's quicker.

 *

 * VIPS uses the fftw Fourier Transform library. If this library was not

 * available when VIPS was configured, these functions will fail.

 *

 * See also: vips_fwfft().

 *

 * Returns: 0 on success, -1 on error.

 */

int

vips_invfft(VipsImage *in, VipsImage **out, ...)

{

  va_list ap;

  int result;

  va_start(ap, out);

  result = vips_call_split("invfft", ap, in, out);

  va_end(ap);

  return result;

}