Mirroring Quaternions

Question

For testing my quaternion averaging algorithm I would like to setup a function that generates an amount of random quaternions where the mean is known. Feeding these quaternions to my function should lead then to the same result as expected by the testing function.

My idea was to represent the quaternion in the angle axis representation (source www.eater.net/quaternions/video/intro).

$q=\cos(\theta)+\sin(\theta)(ai+bj+ck)$

With this representation I thought to generate a random normalized quaternion by getting three random numbers for $\theta$, $a$, $b$ and $c$. then just rotate it about an angle of $(\theta+\pi/2)$. Does this make sense? My quaternions do not turn out the right way. Could you give me a hint where my misunderstanding is?

Here I also got the function in C++ where I want to mirror my quaternions around the quaternion(1,0,0,0) defined in WXYZ:

void getTwoMirroredQuaternions(Eigen::Vector4d& quaternion1,
                                         Eigen::Vector4d& quaternion2) {
  double PI = 3.14159265359;

  double i = getRandomNumber();
  double j = getRandomNumber();
  double k = getRandomNumber();

  i = i/(pow(i,2) + pow(j,2) + pow(k,2));
  j = j/(pow(i,2) + pow(j,2) + pow(k,2));
  k = k/(pow(i,2) + pow(j,2) + pow(k,2));

  std::cout << "i: " << i << " j: " << j << " k: " << k << std::endl;

  double realAngle = getRandomNumber();

  double realPart1 = std::cos(realAngle);
  double iPart1 = std::sin(realAngle)*i;
  double jPart1 = std::sin(realAngle)*j;
  double kPart1 = std::sin(realAngle)*k;

  double realPart2 = std::cos(realAngle+PI/2);
  double iPart2 = std::sin(realAngle+PI/2)*i;
  double jPart2 = std::sin(realAngle+PI/2)*j;
  double kPart2 = std::sin(realAngle+PI/2)*k;

  quaternion1.x() = iPart1;
  quaternion1.y() = jPart1;
  quaternion1.z() = kPart1;
  quaternion1.w() = realPart1;

  quaternion2.x() = iPart2;
  quaternion2.y() = jPart2;
  quaternion2.z() = kPart2;
  quaternion2.w() = realPart2;

  quaternion1.normalize();
  quaternion2.normalize();
}