https://doi.org/10.5201/ipol.2016.172
Raw File
Tip revision: edce3fbb7ce8b72fd316f77f13ff3f1d75091a52 authored by Software Heritage on 19 February 2016, 00:00:00 UTC
ipol: Deposit 1282 in collection ipol
Tip revision: edce3fb
bicubic_interpolation.h
// This program is free software: you can use, modify and/or redistribute it
// under the terms of the simplified BSD License. You should have received a
// copy of this license along this program. If not, see
// <http://www.opensource.org/licenses/bsd-license.html>.
//
// Copyright (C) 2012, Javier Sánchez Pérez <jsanchez@dis.ulpgc.es>
// Copyright (C) 2014, Nelson Monzón López <nmonzon@ctim.es>
// All rights reserved.

#ifndef BICUBIC_INTERPOLATION_H
#define BICUBIC_INTERPOLATION_H

#define BOUNDARY_CONDITION 0

//0 Neumann
//1 Periodic
//2 Symmetric

/**
  *
  * Neumann boundary condition test
  *
**/
int
neumann_bc (int x, int nx, bool & out)
{
  if (x < 0)
    {
      x = 0;
      out = true;
    }
  else if (x >= nx)
    {
      x = nx - 1;
      out = true;
    }

  return x;
}

/**
  *
  * Periodic boundary condition test
  *
**/
int
periodic_bc (int x, int nx, bool & out)
{
  if (x < 0)
    {
      const int n = 1 - (int) (x / (nx + 1));
      const int ixx = x + n * nx;

      x = ixx % nx;
      out = true;
    }
  else if (x >= nx)
    {
      x = x % nx;
      out = true;
    }

  return x;
}


/**
  *
  * Symmetric boundary condition test
  *
**/
int
symmetric_bc (int x, int nx, bool & out)
{
  if (x < 0)
    {
      const int border = nx - 1;
      const int xx = -x;
      const int n = (int) (xx / border) % 2;

      if (n)
	x = border - (xx % border);
      else
	x = xx % border;
      out = true;
    }
  else if (x >= nx)
    {
      const int border = nx - 1;
      const int n = (int) (x / border) % 2;

      if (n)
	x = border - (x % border);
      else
	x = x % border;
      out = true;
    }

  return x;
}


/**
  *
  * Bicubic interpolation in one dimension
  *
**/
double
cubic_interpolation (double v[4],	//interpolation points
		     double x		//point to be interpolated
  )
{
  return v[1] + 0.5 * x * (v[2] - v[0]
			   + x * (2.0 * v[0] - 5.0 * v[1] + 4.0 * v[2] - v[3]
				  + x * (3.0 * (v[1] - v[2]) + v[3] - v[0])));
}


/**
  *
  * Bicubic interpolation in two dimension
  *
**/
double
bicubic_interpolation (double p[4][4],	//array containing the interpolation points
		       double x,	//x position to be interpolated
		       double y		//y position to be interpolated
  )
{
  double v[4];
  v[0] = cubic_interpolation (p[0], y);
  v[1] = cubic_interpolation (p[1], y);
  v[2] = cubic_interpolation (p[2], y);
  v[3] = cubic_interpolation (p[3], y);
  return cubic_interpolation (v, x);
}

/**
  *
  * Compute the bicubic interpolation of a point in an image. 
  * Detects if the point goes outside the image domain
  *
**/
float
bicubic_interpolation (const float *input,	//image to be interpolated
		       const float uu,		//x component of the vector field
		       const float vv,		//y component of the vector field
		       const int nx,		//width of the image
		       const int ny,		//height of the image
		       const int nz,            //number of channels of the image
		       const int k,  		//actual channel
		       const bool border_out = false	//if true, put zeros outside the region
  )
{
  const int sx = (uu < 0) ? -1 : 1;
  const int sy = (vv < 0) ? -1 : 1;

  int x, y, mx, my, dx, dy, ddx, ddy;
  bool out = false;

  switch (BOUNDARY_CONDITION)
    {

    case 0:
      x = neumann_bc ((int) uu, nx, out);
      y = neumann_bc ((int) vv, ny, out);
      mx = neumann_bc ((int) uu - sx, nx, out);
      my = neumann_bc ((int) vv - sx, ny, out);
      dx = neumann_bc ((int) uu + sx, nx, out);
      dy = neumann_bc ((int) vv + sy, ny, out);
      ddx = neumann_bc ((int) uu + 2 * sx, nx, out);
      ddy = neumann_bc ((int) vv + 2 * sy, ny, out);
      break;

    case 1:
      x = periodic_bc ((int) uu, nx, out);
      y = periodic_bc ((int) vv, ny, out);
      mx = periodic_bc ((int) uu - sx, nx, out);
      my = periodic_bc ((int) vv - sx, ny, out);
      dx = periodic_bc ((int) uu + sx, nx, out);
      dy = periodic_bc ((int) vv + sy, ny, out);
      ddx = periodic_bc ((int) uu + 2 * sx, nx, out);
      ddy = periodic_bc ((int) vv + 2 * sy, ny, out);
      break;

    case 2:
      x = symmetric_bc ((int) uu, nx, out);
      y = symmetric_bc ((int) vv, ny, out);
      mx = symmetric_bc ((int) uu - sx, nx, out);
      my = symmetric_bc ((int) vv - sx, ny, out);
      dx = symmetric_bc ((int) uu + sx, nx, out);
      dy = symmetric_bc ((int) vv + sy, ny, out);
      ddx = symmetric_bc ((int) uu + 2 * sx, nx, out);
      ddy = symmetric_bc ((int) vv + 2 * sy, ny, out);
      break;

    default:
      x = neumann_bc ((int) uu, nx, out);
      y = neumann_bc ((int) vv, ny, out);
      mx = neumann_bc ((int) uu - sx, nx, out);
      my = neumann_bc ((int) vv - sx, ny, out);
      dx = neumann_bc ((int) uu + sx, nx, out);
      dy = neumann_bc ((int) vv + sy, ny, out);
      ddx = neumann_bc ((int) uu + 2 * sx, nx, out);
      ddy = neumann_bc ((int) vv + 2 * sy, ny, out);
      break;
    }

  if (out && border_out)

    return 0.0;

  else
    {
      //obtain the interpolation points of the image
      const float p11 = input[(mx  + nx * my) * nz + k];
      const float p12 = input[(x   + nx * my) * nz + k];
      const float p13 = input[(dx  + nx * my) * nz + k];
      const float p14 = input[(ddx + nx * my) * nz + k];

      const float p21 = input[(mx  + nx * y) * nz + k];
      const float p22 = input[(x   + nx * y) * nz + k];
      const float p23 = input[(dx  + nx * y) * nz + k];
      const float p24 = input[(ddx + nx * y) * nz + k];

      const float p31 = input[(mx  + nx * dy) * nz + k];
      const float p32 = input[(x   + nx * dy) * nz + k];
      const float p33 = input[(dx  + nx * dy) * nz + k];
      const float p34 = input[(ddx + nx * dy) * nz + k];

      const float p41 = input[(mx  + nx * ddy) * nz + k];
      const float p42 = input[(x   + nx * ddy) * nz + k];
      const float p43 = input[(dx  + nx * ddy) * nz + k];
      const float p44 = input[(ddx + nx * ddy) * nz + k];

      //create array
      double pol[4][4] = { 
	{p11, p21, p31, p41}, {p12, p22, p32, p42},
	{p13, p23, p33, p43}, {p14, p24, p34, p44}
      };

      //return interpolation
      return bicubic_interpolation (pol, (float) uu - x, (float) vv - y);
    }
}


/**
  *
  * Compute the bicubic interpolation of an image.
  *
**/
void
bicubic_interpolation (const float *input,	//image to be warped
		       const float *u,		//x component of the vector field
		       const float *v,		//y component of the vector field
		       float *output,		//warped output image with bicubic interpolation
		       const int nx,		//width of the image
		       const int ny,		//height of the image
		       const int nz,		//number of channels of the image	
		       bool border_out = false	//if true, put zeros outside the region
  )
{
  for(int k = 0; k < nz; k++){
    #pragma omp parallel for
    for (int i = 0; i < ny; i++)
       for (int j = 0; j < nx; j++){

	      const int p = i * nx + j;
	      const float uu = (float) (j + u[p]);
	      const float vv = (float) (i + v[p]);

	      //obtain the bicubic interpolation at position (uu, vv)
	      output[p * nz + k] = bicubic_interpolation (input, uu, vv, nx, ny, nz, k, border_out);
      
      }
  } // end multi-channel loop
}

#endif
back to top