Why my cosine interpolation of a cube's face doesn't work?

Question

NB please : executable use cases are available at the end of this question. I begin this question by showing you the problems of my program, then I explain how the latter works, and finally I end up with these both executable use cases.

What I want to do

I want to draw the bottom face of a cube. To do that, I interpolate between the 4 vertices of this face. I use two kinds of interpolation :

1. Linear

2. Cosine

I know it's typically the aim of GPUs, but I wanted to implement it myself to learn and try the "deep" mechanisms of 3D rendering.

Results

Each point of the following both pictures are interpolated between the 4 vertices of the bottom face (except the 8 vertices of my cube).

With linear interpolation

It's OK, no problem.

With cosine interpolation

Error. As you can see, there is some weird problem : face's points are inclined towards left.

How do I interpolate between 4 points ?

I test all the combinations of interpolation's weights and, if the sum of the weights equals 1, I draw the resulting point. I know this method is not very good (a lot of consumption of CPU and RAM), and there are P duplicates. But this problem is (1) secondary and (2) out of the range of this StackOverflow question.

In other words : the definition of the interpolation between 4 points is : P = aA + bB + cC + dD. a + b + c + d must equals 1 to draw P, otherwise it's not drawn. (a;b;c;d) \in [0;1]^4 by the way.

A concrete example is (step = 1) :

P = 0a + 0b + 0c + 0d is NOT drawn

P = 0a + 0b + 0c + 1d is drawn

P = 0a + 0b + 1c + 0d is drawn

P = 0a + 0b + 1c + 1d is NOT drawn

etc. ("I test all the combinations"). (Last is : P = 1a + 1b + 1c + 1d, which is not drawn).

To test all these combinations, I iterates on the weights, with a step, using recursion.

Case of the linear interpolation

The above definition doesn't change : P = aA + bB + cC + dD.

Case of the cosine interpolation

The above definition does change : P = aA + bB + cC + dD becomes P = a'A + b'B + c'C + d'D, with : a' = 0.5(1 - cos(a * \pi)) and the idea is the same for the three other weights.

This new definition is not random, I could explain how I found it, but it's out of the range of this topic (in résumé : a simple remap of the cosine function).

Question

Why does it works well with the linear interpolation, and not the cosine one ? Indeed, the cosine face's points are inclined towards left.

By the way : interpolations between two points work well

Linear and cosine interpolations both work very well when used between two points, as you can see below.

Linear interpolation

It's OK, no problem.

Cosine interpolation

It's OK, no problem.

Implementation and Executable

I wrote the program in Scala. First I show you the functions, then two use cases you can execute along with the given functions.

Tests all the combinations of interpolations between n points of k coordinates

def computeAllPossibleInterpolatedPoints(step : Double, points : Seq[Seq[Double]], transform: (Double) => Double) : Seq[Seq[Double]] = {
  var returned_object : Seq[Seq[Double]] = Seq.empty[Seq[Double]]
  recursiveInterpolation(0, Seq.empty[Double])

  def recursiveInterpolation(current_weight_id : Int, building_weights : Seq[Double]) : Unit = {

    (.0 to 1.0 by step).foreach(current_step => {
      if (current_weight_id < points.size - 1) {
        recursiveInterpolation(current_weight_id + 1, building_weights :+ current_step)
      } else {
        val found_solution = (building_weights :+ current_step).map(transform)
        if(BigDecimal(found_solution.sum).setScale(5, BigDecimal.RoundingMode.HALF_UP).toDouble == 1.0) {
          returned_object = returned_object :+ interpolation(found_solution, points)
        }
      }
    })
  }

  returned_object
}

Interpolates between n points of k coordinates

def interpolation(weights: Seq[Double], points: Seq[Seq[Double]], transform: (Double) => Double = null) : Seq[Double] = {
  var transformed_weights = weights
  if(transform != null) {
    transformed_weights = weights.map(transform)
  }
  if(BigDecimal(transformed_weights.sum).setScale(5, BigDecimal.RoundingMode.HALF_UP).toDouble != 1.0) {
    println("ERROR : `SUM(weights) != 1`. Returning `null`.")
    return null
  }
  if(transformed_weights.exists(weight => BigDecimal(weight).setScale(5, BigDecimal.RoundingMode.HALF_UP).toDouble < 0)
    ||
    transformed_weights.exists(weight => BigDecimal(weight).setScale(5, BigDecimal.RoundingMode.HALF_UP).toDouble > 1)) {
    println("ERROR : `EXISTS(weight) / weight < -1 OR weight > 1`. Returning `null`.")
    return null
  }

  transformed_weights.zip(points).map(
    weight_point => weight_point._2.map(coordinate => weight_point._1 * coordinate)
  ).reduce((point_a : Seq[Double], point_b : Seq[Double]) => point_a.zip(point_b).map(coordinate_points => coordinate_points._1 + coordinate_points._2))
}

Interpolation functions

def linear(weight : Double) : Double = {
  weight
}

def cosine(weight : Double) : Double = {
  (1 - Math.cos(weight * Math.PI)) * 0.5
}

Two use cases : between 2 points, each made of 2 coordinates (first use case : linear interpolation, second use case : cosine interpolation)

Note that you can add 2 new points to each of the following functions, to interpolate a face instead of a segment (for now, the below code indeed interpolates points of a segment).

val interpolated_points : Seq[Seq[Int]] = interpolator.computeAllPossibleInterpolatedPoints(0.05, Seq(Seq(22, 22), Seq(33, 33)), interpolator.linear).map(_.map(_.intValue()))

val interpolated_points : Seq[Seq[Int]] = interpolator.computeAllPossibleInterpolatedPoints(0.05, Seq(Seq(22, 22), Seq(33, 33)), interpolator.cosine).map(_.map(_.intValue()))

Answer 1

I only remember after struggling to read through code, but I'd recommend posting pseudo-code rather than code here. Regardless, I tried to understand your implementation, and as far as I can tell, you test weights $a,b,c,d$ such that they sum to 1, and for these $a+b+c+d=1$ you directly plug them into your cosine interpolation to compute $a',b',c',d'$. In general, this will not give you a point in the face because you actually want the constraint $a'+b'+c'+d'=1$ to perform an interpolation.

Now, why does this work with a segment interpolation and not a face interpolation? The tl;dr is that $a+b=1$ will ensure $a'+b'=1$. Indeed $$\cos(a\pi)=\cos(\pi-b\pi)=-\cos(-b\pi)=-\cos(b\pi)$$ and this yields $$a'+b' = \frac 12(1-\cos(a\pi)) +\frac 12(1-\cos(b\pi))=1$$ When you have more coefficients however, this is unlikely to hold, so your interpolation is not really an interpolation.

---

To solve this, you can use the "transform" (linear or cosine functions) before testing the validity of the weights.