int ViewElemInterpolate()

in layer1/View.cpp [805:1291]


int ViewElemInterpolate(PyMOLGlobals * G, CViewElem * first, CViewElem * last,
                        float power, float bias,
                        int simple, float linearity, int hand, float cut)
{
  float first3x3[9];
  float last3x3[9];
  float inverse3x3[9];
  float inter3x3[9];
  float rot_axis[3], trans_axis[3] = { 0.0F, 0.0F, 0.0F };
  float angle;
  CViewElem *current;
  ov_diff n = (last - first) - 1;
  Matrix53f rot, imat;
  int a;
  float tVector[3], tCenter[3], tDir[3];
  float tLen = 0.0F;
  float bisect[3], v2[3];
  float translate_angle = 0.0F;
  float pivot[3] = { 0.0F, 0.0F, 0.0F };
  const float _1 = 1.0F, _p5 = 0.5F;
  int parabolic = true;
  int timing_flag;
  double timing = 0.0F;
  int state_flag;
  int state = 0;
  float pre[3];
  float firstC44f[16], firstRTTT[16], firstR44f[16];
  float lastC44f[16], lastRTTT[16], lastR44f[16];
  int linear = false;
  int debug = Feedback(G,FB_Movie, FB_Debugging);
  
  if(hand == 0)
    hand = 1;

  if(debug) {
    printf("ViewElemInterpolate: %8.3f %8.3f %d %8.3f %d %8.3f\n",
           power, bias, simple, linearity, hand, cut);
    dump44d(first->matrix,"first->matrix");
    dump44d(last->matrix,"last->matrix");
    printf("first->pre_flag %d first->post_flag %d\n",first->pre_flag, first->post_flag);
    dump3d(first->pre,"first->pre");
    dump3d(first->post,"first->post");
    printf("last->pre_flag %d last->post_flag %d\n",last->pre_flag, last->post_flag);
    dump3d(last->pre,"last->pre");
    dump3d(last->post,"last->post");
  }
  if(power == 0.0F) {
    if(first->power_flag && last->power_flag) {
      if(((first->power > 0.0F) && (last->power > 0.0F)) ||
         ((first->power < 0.0F) && (last->power < 0.0F))) {
        power = (first->power + last->power) / 2.0F;
      } else if(fabs(first->power) > fabs(last->power)) {
        power = first->power;
      } else if(last->power < 0.0F) {
        power = last->power;
      } else {
        power = first->power;
      }
    } else if(first->power_flag) {
      power = first->power;
    } else if(last->power_flag) {
      power = last->power;
    } else {
      power = 1.4F; /* default */
    }
  }
  if(power < 0.0F) {
    parabolic = false;
    power = -power;
  }

  if(bias < 0.0F) { /* default */
    if(first->bias_flag && last->bias_flag) {
      if((first->bias > 0.0F) && (last->bias > 0.0F)) {
        bias = (first->bias * 1.0F/last->bias);
      } else if(fabs(first->bias) > 0.0) {
        bias = first->bias;
      } else if(last->bias > 0.0F) {
        bias = 1.0F/last->bias;
      } else {
        bias = 1.0F;
      }
    } else if(first->bias_flag) {
      bias = first->bias;
    } else if(last->bias_flag) {
      bias = 1.0F/last->bias;
    } else {
      bias = 1.0F; /* default */
    }
  }

  if(bias <= 0.0F) {
    bias = 1.0F;
  }

  /* WARNING: this routine is operating on column-major matrices!!! */

  copy44d33f(first->matrix, first3x3);
  copy44d33f(last->matrix, last3x3);

  transpose33f33f(first3x3, inverse3x3);

  multiply33f33f(inverse3x3, last3x3, &rot[0][0]);      /* [rot] = [first]^-1 [last] */
  matrix_to_rotation(rot, rot_axis, &angle);

  if(debug)
    dump3f(rot_axis, "rot_axis");

  if(hand) {
    if((cPI - fabs(angle)) < 0.01F) {   /* this a complete 180 degree motion */
      if(((rot_axis[0] * 0.7F + rot_axis[1] * 0.8F + rot_axis[2] * 0.9F) * hand * angle) >
         0.0F) {
        invert3f(rot_axis);
        if(angle > 0) {
          angle = (float) ((2 * cPI) - angle);
        } else {
          angle = (float) (-(2 * cPI) - angle);
        }
      }
    }
  }

  if(!simple) {
    /* switch back into row major to promote developer sanity */

    copy33f44f(first3x3, firstC44f);
    copy33f44f(last3x3, lastC44f);

    transpose44f44f(firstC44f, firstRTTT);
    transpose44f44f(lastC44f, lastRTTT);

    /* form TTTs */

    firstRTTT[12] = (float) -first->pre[0];
    firstRTTT[13] = (float) -first->pre[1];
    firstRTTT[14] = (float) -first->pre[2];

    firstRTTT[3] = (float) first->post[0];
    firstRTTT[7] = (float) first->post[1];
    firstRTTT[11] = (float) first->post[2];

    lastRTTT[12] = (float) -last->pre[0];
    lastRTTT[13] = (float) -last->pre[1];
    lastRTTT[14] = (float) -last->pre[2];

    lastRTTT[3] = (float) last->post[0];
    lastRTTT[7] = (float) last->post[1];
    lastRTTT[11] = (float) last->post[2];

    if(debug)
      dump44f(firstRTTT, "firstRTTT");
    if(debug)
      dump44f(lastRTTT, "lastRTTT");

    /* convert to homogenous */

    convertTTTfR44f(firstRTTT, firstR44f);
    convertTTTfR44f(lastRTTT, lastR44f);

    /* reset both matrices to a common origin */

    {
      float first_pre[3], last_pre[3];
      float post[4], *dst;

      copy3d3f(first->pre, first_pre);
      copy3d3f(last->pre, last_pre);
      average3f(first_pre, last_pre, pre);

      transform44f3fas33f3f(firstR44f, pre, post);
      copy44f(firstR44f, firstRTTT);
      firstRTTT[3] += post[0];
      firstRTTT[7] += post[1];
      firstRTTT[11] += post[2];
      dst = firstRTTT + 12;
      invert3f3f(pre, dst);

      transform44f3fas33f3f(lastR44f, pre, post);
      copy44f(lastR44f, lastRTTT);
      lastRTTT[3] += post[0];
      lastRTTT[7] += post[1];
      lastRTTT[11] += post[2];
      dst = lastRTTT + 12;
      invert3f3f(pre, dst);
    }

    if(debug)
      dump44f(firstRTTT, "firstRTTT");
    if(debug)
      dump44f(lastRTTT, "lastRTTT");

    /*    convertTTTfR44f(firstRTTT, firstR44f);
       convertTTTfR44f(lastRTTT, lastR44f); */

    /* now populate the translation fields */

    rot[3][0] = firstRTTT[3];
    rot[3][1] = firstRTTT[7];
    rot[3][2] = firstRTTT[11];

    rot[4][0] = lastRTTT[3];
    rot[4][1] = lastRTTT[7];
    rot[4][2] = lastRTTT[11];

    /* now set up the interpolation */

    subtract3f(&rot[4][0], &rot[3][0], tVector);
    tLen = (float) length3f(tVector);
    average3f(&rot[4][0], &rot[3][0], tCenter);

    if(tLen < 0.0001F) {
      if(debug)
        printf("translation too short %8.3f\n", tLen);
      simple = true;
    }
  }

  if(!simple) {

    normalize23f(tVector, tDir);
    if(debug)
      dump3f(tDir, "tDir");
    cross_product3f(tDir, rot_axis, bisect);
    /* bisect is a vector in the translation arc */
    if(length3f(bisect) < 0.0001F) {
      if(debug)
        printf("rotation coincident with translation\n");
      linear = true;
    }
  }

  if(!(simple || linear)) {
    normalize3f(bisect);

    /* this section needs work... */

    cross_product3f(bisect, tDir, trans_axis);
    normalize3f(trans_axis);

    transform33Tf3f(&rot[0][0], bisect, v2);    /* column major */

    remove_component3f(v2, trans_axis, v2);
    normalize3f(v2);            /* project vector onto plane _|_ to axis */

    if(debug) {
      dump3f(rot_axis, "rot_axis");
      dump3f(tDir, "tDir");
      dump3f(bisect, "bisect");
      dump3f(trans_axis, "trans_axis");
      dump3f(v2, "v2");
    }

    {
      double dot = dot_product3f(bisect, v2);
      if(dot < -1.0F)
        dot = -1.0F;
      if(dot > 1.0F)
        dot = 1.0F;
      translate_angle = (float) acos(dot);

      /* if translation angle > rotation angle then sets translation angle 
       * to same as rotation angle, with proper sign of course */

      if((fabs(translate_angle) > fabs(angle)) && (fabs(angle) > R_SMALL4)) {
        translate_angle = (float) (fabs(angle) * (translate_angle / fabs(angle)));
      }

      if(fabs(translate_angle) < 0.0001F) {
        linear = true;
        if(debug)
          printf("no significant rotation\n");
      }

      if((translate_angle * angle) < 0.0F) {
        /* if motions are in opposing directions, then flip translation axis and location */
        invert3f(bisect);
        invert3f(trans_axis);
      }

    }
  }

  if(!(simple || linear)) {
    float pLen = (float) tan(translate_angle / 2);
    if(fabs(pLen) > 0.0000001)
      pLen = (tLen / 2) / pLen;
    else {
      if(debug)
        printf("pLen too short %8.3f\n", pLen);
      simple = true;
    }

    if(!simple) {
      pivot[0] = tCenter[0] + pLen * bisect[0];
      pivot[1] = tCenter[1] + pLen * bisect[1];
      pivot[2] = tCenter[2] + pLen * bisect[2];
    }

    if(debug && !simple) {
      dump3f(tCenter, "center");
      dump3f(pivot, "pivot");
      printf("pLen %8.3f angle %8.3f translate_angle %8.3f\n",
             pLen, angle, translate_angle);
    }
  }

  /* now interpolate */

  state_flag = first->state_flag && last->state_flag;

  timing_flag = first->timing_flag && last->timing_flag;

  current = first + 1;

  if(debug)
    dump44f(firstR44f, "first");

  for(a = 0; a < n; a++) {
    float fxn = ((float) a + 1) / (n + 1);
    float fxn_1;

    if(timing_flag) {
      timing = (first->timing * (1.0F - fxn) + (last->timing * fxn));
    }

    if(state_flag) { /* states are interpolated linearly by default */
      state = (int)(first->state * (1.0F - fxn) + (last->state * fxn) + 0.499F);
    }

    if(bias != 1.0F) {
      fxn = 1 - (float) pow(1 - pow(fxn, bias), _1 / bias);
    }

    if((power != 1.0F) || (!parabolic)) {
      if(fxn < 0.5F) {
        if(!parabolic)
          fxn = (float) ((_1 - cos(cPI * fxn)) * _p5);  /* circular */
        fxn = (float) pow(fxn * 2.0F, power) * _p5;     /* parabolic */
      } else if(fxn > 0.5F) {
        fxn = _1 - fxn;
        if(!parabolic)
          fxn = (float) ((_1 - cos(cPI * fxn)) * _p5);
        fxn = (float) pow(fxn * 2.0F, power) * _p5;     /* parabolic */
        fxn = _1 - fxn;
      }
    }

    fxn_1 = 1.0F - fxn;

    ViewElemCopy(G, first, current);

    if(simple) {
      rotation_matrix3f(fxn * angle, rot_axis[0], rot_axis[1], rot_axis[2], &imat[0][0]);


/* [cur] = [first] [partial-rot], so....
   at start: [cur] = [first] [identity] = [first]
   at end: [cur] = [first] [first]^-1 [last] = [last] 
*/
      current->matrix_flag = true;
      multiply33f33f(first3x3, &imat[0][0], inter3x3);

      copy33f44d(inter3x3, current->matrix);

      if(first->pre_flag && last->pre_flag) {
        mix3d(first->pre, last->pre, (double) fxn, current->pre);
        current->pre_flag = true;
      } else {
        current->pre_flag = false;
      }
      if(first->post_flag && last->post_flag) {
        mix3d(first->post, last->post, (double) fxn, current->post);
        current->post_flag = true;
      } else {
        current->post_flag = false;
      }
    } else if(linear) {
      int b;
      rotation_matrix3f(fxn * angle, rot_axis[0], rot_axis[1], rot_axis[2], &imat[0][0]);
      current->matrix_flag = true;
      multiply33f33f(first3x3, &imat[0][0], inter3x3);

      copy33f44d(inter3x3, current->matrix);

      current->pre_flag = true;
      copy3f3d(pre, current->pre);

      current->post_flag = true;
      for(b = 0; b < 3; b++) {
        imat[4][b] = (float) ((1.0 - fxn) * rot[3][b] + fxn * rot[4][b]);
      }
      copy3f3d(&imat[4][0], current->post);
    } else {
      matrix_interpolate(imat, rot,
                         pivot, bisect,
                         rot_axis, angle, trans_axis, translate_angle, fxn, linearity);

      current->matrix_flag = true;
      multiply33f33f(first3x3, &imat[0][0], inter3x3);

      copy33f44d(inter3x3, current->matrix);

      current->pre_flag = true;
      copy3f3d(pre, current->pre);

      current->post_flag = true;
      copy3f3d(&imat[4][0], current->post);

    }
    if(debug) {
      if((a == 0) || (a == n - 1)) {
        float curC44f[16], curRTTT[16], curR44f[16];

        copy33f44f(inter3x3, curC44f);

        transpose44f44f(curC44f, curRTTT);

        /* form TTTs */

        curRTTT[12] = (float) -current->pre[0];
        curRTTT[13] = (float) -current->pre[1];
        curRTTT[14] = (float) -current->pre[2];

        curRTTT[3] = (float) current->post[0];
        curRTTT[7] = (float) current->post[1];
        curRTTT[11] = (float) current->post[2];

        convertTTTfR44f(curRTTT, curR44f);
        dump44f(curR44f, "cur");
      }
    }

    if(first->clip_flag && last->clip_flag) {
      current->front = first->front * fxn_1 + last->front * fxn;
      current->back = first->back * fxn_1 + last->back * fxn;
      current->clip_flag = true;
    } else {
      current->clip_flag = false;
    }

    if(first->ortho_flag && last->ortho_flag) {
      float approx_ortho = first->ortho * fxn_1 + last->ortho * fxn;
      if(first->pre_flag && last->pre_flag) {
        float first_far = first->pre[2] * tan(cPI * fabs(first->ortho) / 360.0);
        float last_far = last->pre[2] * tan(cPI * fabs(last->ortho) / 360.0);

        float cur_far = first_far * fxn_1 + last_far * fxn;
        current->ortho = 360.0 * atan(cur_far / current->pre[2]) / cPI;

        if((current->ortho * approx_ortho) < 0) /* fix sign */
          current->ortho = -current->ortho;
      } else {
        current->ortho = approx_ortho;
      }
    }
    current->specification_level = 1;

    if(state_flag) {
      current->state_flag = true;
      current->state = state;
    }

    if(timing_flag) {
      current->timing_flag = true;
      current->timing = timing;
    }


    if(first->scene_flag && last->scene_flag) {
      if(current->scene_name) {
        OVLexicon_DecRef(G->Lexicon, current->scene_name);
      }
      current->scene_flag = true;
      if(fxn >= cut) {
        current->scene_name = last->scene_name;
      } else {
        current->scene_name = first->scene_name;
      }
      OVLexicon_IncRef(G->Lexicon, current->scene_name);
    }
    current++;
  }
  if(debug)
    dump44f(lastR44f, "last");

  return 1;
}