/* forward FFT

 *

 * Author: Nicos Dessipris

 * Written on: 12/04/1990

 * Modified on : 09/05/1990 to cope with float input

 * Modified on : 08/03/1991   history removed

 * Modified on : 03/04/1991 to cope with any input

 *

 * 28/6/95 JC

 *  - rewritten to use im_clip2f() rather than own code

 *  - memory leaks fixed

 * 10/9/98 JC

 *  - frees memory more quickly

 * 2/4/02 JC

 *  - fftw code added

 * 13/7/02 JC

 *  - output Type set to IM_TYPE_FOURIER to help nip

 * 27/2/03 JC

 *  - exploits real_to_complex() path in libfftw for real input (thanks

 *    Matt) for a 2x speed-up

 * 17/11/03 JC

 *  - fix a segv for wider than high images in the real_to_complex() path

 *    (thanks Andrey)

 *  - fixes to real_to_complex() path to give the correct result for

 *    non-square images, including odd widths and heights

 * 3/11/04

 *  - added fftw3 support

 * 7/2/10

 *  - cleanups

 *  - gtkdoc

 * 25/3/10

 *  - have a "t" image linked to out to keep the image alive for longer

 * 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 _VipsFwfft {

  VipsFreqfilt parent_instance;

} VipsFwfft;

typedef VipsFreqfiltClass VipsFwfftClass;

G_DEFINE_TYPE(VipsFwfft, vips_fwfft, VIPS_TYPE_FREQFILT);

/* Everything in fftw3 except execute has to be behind a mutex.

 */

GMutex *vips__fft_lock = NULL;

static void *

vips__fft_thread_init(void *data)

{

  vips__fft_lock = vips_g_mutex_new();

    return NULL;

}

void

vips__fft_init(void)

{

  static GOnce once = G_ONCE_INIT;

  VIPS_ONCE(&once, vips__fft_thread_init, NULL);

}

/* Real to complex forward transform.

 */

static int

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

{

  VipsFwfft *fwfft = (VipsFwfft *) object;

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

  VipsObjectClass *class = VIPS_OBJECT_GET_CLASS(fwfft);

  const guint64 size = VIPS_IMAGE_N_PELS(in);

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

  double *half_complex;

  double *planner_scratch;

  fftw_plan plan;

  double *buf, *q, *p;

  int x, y;

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

    vips_check_uncoded(class->nickname, in))

    return -1;

  /* Convert input to a real double membuffer.

   */

  t[1] = vips_image_new_memory();

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

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

    return -1;

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

   * height and ny for width. Use a separate scratch buffer for the

   * planner, we can't overwrite real->data

   */

  if (!(planner_scratch = VIPS_ARRAY(fwfft,

        VIPS_IMAGE_N_PELS(in), double)))

    return -1;

  if (!(half_complex = VIPS_ARRAY(fwfft,

        in->Ysize * half_width * 2, double)))

    return -1;

  g_mutex_lock(vips__fft_lock);

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

        planner_scratch, (fftw_complex *) half_complex,

        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_r2c(plan,

    (double *) t[1]->data, (fftw_complex *) half_complex);

  g_mutex_lock(vips__fft_lock);

  fftw_destroy_plan(plan);

  g_mutex_unlock(vips__fft_lock);

  /* Write to out as another memory buffer.

   */

  *out = vips_image_new_memory();

  if (vips_image_pipelinev(*out, VIPS_DEMAND_STYLE_ANY, in, NULL))

    return -1;

  (*out)->BandFmt = VIPS_FORMAT_DPCOMPLEX;

  (*out)->Type = VIPS_INTERPRETATION_FOURIER;

  if (!(buf = VIPS_ARRAY(fwfft, VIPS_IMAGE_N_PELS(*out), double)))

    return -1;

  /* Copy and normalise. The right half is the up/down and

   * left/right flip of the left, but conjugated. Do the first

   * row separately, then mirror around the centre row.

   */

  p = half_complex;

  q = buf;

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

    q[0] = p[0] / size;

    q[1] = p[1] / size;

    p += 2;

    q += 2;

  }

  p = half_complex + ((in->Xsize + 1) / 2 - 1) * 2;

  for (x = half_width; x < (*out)->Xsize; x++) {

    q[0] = p[0] / size;

    q[1] = -1.0 * p[1] / size;

    p -= 2;

    q += 2;

  }

  if (vips_image_write_line(*out, 0, (VipsPel *) buf))

    return -1;

  for (y = 1; y < (*out)->Ysize; y++) {

    p = half_complex + y * half_width * 2;

    q = buf;

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

      q[0] = p[0] / size;

      q[1] = p[1] / size;

      p += 2;

      q += 2;

    }

    /* Good grief.

     */

    p = half_complex + 2 * /* clang-format off */

        (((*out)->Ysize - y + 1) * half_width - 2 +

          (in->Xsize & 1));

    /* clang-format on */

    for (x = half_width; x < (*out)->Xsize; x++) {

      q[0] = p[0] / size;

      q[1] = -1.0 * p[1] / size;

      p -= 2;

      q += 2;

    }

    if (vips_image_write_line(*out, y, (VipsPel *) buf))

      return -1;

  }

  return 0;

}

/* Complex to complex forward transform.

 */

static int

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

{

  VipsFwfft *fwfft = (VipsFwfft *) object;

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

  VipsObjectClass *class = VIPS_OBJECT_GET_CLASS(fwfft);

  fftw_plan plan;

  double *planner_scratch;

  double *buf, *q, *p;

  int x, y;

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

    vips_check_uncoded(class->nickname, in))

    return -1;

  /* 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;

  /* We have to have a separate buffer for the planner to work on.

   */

  if (!(planner_scratch = VIPS_ARRAY(fwfft,

        VIPS_IMAGE_N_PELS(in) * 2, double)))

    return -1;

  /* Make the plan for the transform.

   */

  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_FORWARD,

        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 *) t[1]->data, (fftw_complex *) t[1]->data);

  g_mutex_lock(vips__fft_lock);

  fftw_destroy_plan(plan);

  g_mutex_unlock(vips__fft_lock);

  /* Write to out as another memory buffer.

   */

  *out = vips_image_new_memory();

  if (vips_image_pipelinev(*out, VIPS_DEMAND_STYLE_ANY, in, NULL))

    return -1;

  (*out)->BandFmt = VIPS_FORMAT_DPCOMPLEX;

  (*out)->Type = VIPS_INTERPRETATION_FOURIER;

  if (!(buf = VIPS_ARRAY(fwfft, VIPS_IMAGE_N_PELS(*out), double)))

    return -1;

  /* Copy to out, normalise.

   */

  p = (double *) t[1]->data;

  for (y = 0; y < (*out)->Ysize; y++) {

    guint64 size = VIPS_IMAGE_N_PELS(*out);

    q = buf;

    for (x = 0; x < (*out)->Xsize; x++) {

      q[0] = p[0] / size;

      q[1] = p[1] / size;

      p += 2;

      q += 2;

    }

    if (vips_image_write_line(*out, y, (VipsPel *) buf))

      return -1;

  }

  return 0;

}

static int

vips_fwfft_build(VipsObject *object)

{

  VipsFreqfilt *freqfilt = VIPS_FREQFILT(object);

  VipsFwfft *fwfft = (VipsFwfft *) object;

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

  VipsImage *in;

  vips__fft_init();

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

    return -1;

  in = freqfilt->in;

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

    return -1;

  in = t[0];

  if (vips_band_format_iscomplex(in->BandFmt)) {

    if (vips__fftproc(VIPS_OBJECT(fwfft), in, &t[1],

        cfwfft1))

      return -1;

  }

  else {

    if (vips__fftproc(VIPS_OBJECT(fwfft), in, &t[1],

        rfwfft1))

      return -1;

  }

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

    return -1;

  return 0;

}

static void

vips_fwfft_class_init(VipsFwfftClass *class)

{

  VipsObjectClass *vobject_class = VIPS_OBJECT_CLASS(class);

  vobject_class->nickname = "fwfft";

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

  vobject_class->build = vips_fwfft_build;

}

static void

vips_fwfft_init(VipsFwfft *fwfft)

{

}

#endif /*HAVE_FFTW*/

/**

 * vips_fwfft: (method)

 * @in: input image

 * @out: (out): output image

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

 *

 * Transform an image to Fourier space.

 *

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

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

 *

 * See also: vips_invfft().

 *

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

 */

int

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

{

  va_list ap;

  int result;

  va_start(ap, out);

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

  va_end(ap);

  return result;

}