/*DDK*************************************************************************/
/*                                                                           */
/* COPYRIGHT (C) Microsoft Corporation, 1989                                 */
/* COPYRIGHT    Copyright (C) 1995 IBM Corporation                           */
/*                                                                           */
/*    The following IBM OS/2 WARP source code is provided to you solely for  */
/*    the purpose of assisting you in your development of OS/2 WARP device   */
/*    drivers. You may use this code in accordance with the IBM License      */
/*    Agreement provided in the IBM Device Driver Source Kit for OS/2. This  */
/*    Copyright statement may not be removed.                                */
/*                                                                           */
/*****************************************************************************/
#pragma  pagesize(55)

/**************************************************************************
 *
 * SOURCE FILE NAME = XFORMS.C
 *
 * DESCRIPTIVE NAME = Transformation Functions
 *
 *
 * VERSION = V2.0
 *
 * DATE      09/18/88
 *
 * DESCRIPTION Conversion/Transfromation functions
 *
 *
 * FUNCTIONS
 *             NotifyTransformChange
 *             get_angle_xform
 *             ConvertWithMatrix
 *             f_lsqrt
 *             newton_f_lsqrt
 *             f_ldivf
 *
 *
 * NOTES
 *
 *
 * STRUCTURES
 *
 * EXTERNAL REFERENCES
 *
 * EXTERNAL FUNCTIONS
 *
*/
#include "plotters.h"
#include "dispatch.h"
#include "error.h"
#include "xforms.h"
#include "prdmath.h"
#include "lockddc.h"

/*
**  Private function definitions
*/

ULONG  f_lsqrt(ULONG);

ULONG  newton_f_lsqrt(ULONG);
ULONG  f_ldivf(ULONG,ULONG);

/***************************************************************************
 *
 * FUNCTION NAME = NotifyTransformChange
 *
 * DESCRIPTION   = Deals with changes to the transformation matrices
 *                 as reported the engine. Returns the result of an
 *                 indirect call to the GreNotifyTransformChange
 *
 * INPUT         = hDC        - device context handle
 *                 ulFlags    - miscellaneous flags
 *                 pDDC       - pointer to device driver context
 *                 ulFunN     - function number
 *
 * OUTPUT        = pXformData - new transform matrix
 *
 * RETURN-NORMAL = NONE
 * RETURN-ERROR  = NONE
 *          
 *                                                               
 *
 **************************************************************************/

LONG  NotifyTransformChange ( HDC hDC, ULONG ulFlags,
                              PNOTIFYTRANSFORMDATA pXformData,
                              PDDC pDDC, ULONG ulFunN )
/*
**  hDC         Handle to DC
**  ulFlags     Miscellaneous flags about the next parameter
**  pXformData  New transform matrix
**  pDDC        Pointer to cookie
**  ulFunN      Operation code
*/

{
  /*
  **        Function is called when the world->device transform is
  **  changed.  Processing is simple:  firstly transform our current
  **  position,  and if it is outside the 16 bit range allowed,
  **  return an error.   If position is OK,  call the engine's
  **  entry point so that it may perform its housekeeping.
  */

  register PPOINTL pp;
  POINTL PtlTemp[2];                   /* For convert to operate on         */
  LONG lResult;
  RECTL RTemp;
  ULONG ulRet;

  if(!(EnterDriver(pDDC)))
    return(GPI_ERROR);

  /*
  **          In case of MEMORYDC Dispatch to Graphics Engine
  */
  if (pDDC->usDCType == OD_MEMORY)
  {
    if ((ulRet = InnerGreSetModelXform (pDDC->hdcMemory, &pXformData->xform,
                                        SX_OVERWRITE,
                                    (ulFunN & 0xFFFF0000) | NGreSetModelXform)))
      ulRet = (*daNotifyTransformChange) (hDC, ulFlags, pXformData, pDDC,
                                            ulFunN);
    LeaveDriver( pDDC );
    return( ulRet );
  }

  pp = &PtlTemp[0];
  *pp = pDDC->DCState.WorldPosition;   /* Where we are now                  */
  *(pp+1) = pDDC->DCState.abnd.ptlRefPoint;

  if (!((*daConvert)(hDC, CVTC_WORLD, CVTC_DEVICE, pp, 2L, pDDC,
                     NGreConvert)))
  {
    /*
    ** No good - proceed no further
    */
    LeaveDriver(pDDC);
    return  GPI_ERROR;
  }

  /*
  **    Verify that the new current position is still within the
  **  device's coordinate space.  The value should be within the
  **  limits of a signed 16 bit field,  unless the plotter is
  **  larger than that.  Hence the following code sets up the limits.
  **  Note that the scale values used mey not be correct - there is
  **  interaction at the plotter between the P1P2 values and the
  **  scale values.
  */
  RTemp.xLeft   = pDDC->pPDevice->CurrentScale.xLeft > -32768 ?
                     -32768 : pDDC->pPDevice->CurrentScale.xLeft;
  RTemp.xRight  = pDDC->pPDevice->CurrentScale.xRight < 32767L ?
                      32767L : pDDC->pPDevice->CurrentScale.xRight;
  RTemp.yBottom = pDDC->pPDevice->CurrentScale.yBottom > -32768 ?
                     -32768 : pDDC->pPDevice->CurrentScale.yBottom;
  RTemp.yTop    = pDDC->pPDevice->CurrentScale.yTop < 32767L ?
                      32767L : pDDC->pPDevice->CurrentScale.yTop;

  if (pp->x < RTemp.xLeft || pp->x > RTemp.xRight || pp->y < RTemp.yBottom ||
     pp->y > RTemp.yTop)
  {
    GplErrSetError(PMERR_COORDINATE_OVERFLOW);
    LeaveDriver(pDDC);
    return  GPI_ERROR;
  }

  /*
  **        Call the next in chain for whatever that wants to do,  and
  **  if all is well,  update our position and return OK.
  */
  lResult = (*daNotifyTransformChange)(hDC, ulFlags, pXformData, pDDC,
                                       ulFunN);

  /*  If success and any actual change in the transformation matrix
  **  attempt to isolate translation and/or scaling values for later use.
  */
  if (lResult && (pXformData->xform.fxM11 != pDDC->DCState.xfLastMatrix.fxM11
               || pXformData->xform.fxM12 != pDDC->DCState.xfLastMatrix.fxM12 ||
                  pXformData->xform.fxM21 != pDDC->DCState.xfLastMatrix.fxM21 ||
                  pXformData->xform.fxM22 != pDDC->DCState.xfLastMatrix.fxM22))
  {

    /*
    **  The pXformData structure contains four point that are
    **  of interest as far as rotation is concerned.
    **
    **       fxM11, fxM12, fxM21, fxM22
    **
    **  Using a line of known initial properties, transform it
    **  according to the information in pXformData.  Then, using
    **  the new values (post-transformation) "deduce" the current
    **  rotation.
    **
    **    if fxM12 NOT zero OR (fxM12 EQU zero AND fxM21 EQU zero) then
    **       rotation is present
    **    else
    **       no rotation is present.
    **
    **    if ROTATION then
    **       transform unit line according to XformData
    **       use new points and XformData to determine angle
    **       of rotation, according to above table
    **
    **  dlr, mgx-sysdev : 25 Apr 1991
    **  *****************************
    **  This section had special-cased rotation in all 4 quadrants and
    **    their boundaries.  It was producing incorrect results, mirrored
    **    across the x-axis.  This was mostly hidden by errors found in
    **    get_angle_xform.  The new implementation is identical to that
    **    found in the 1.2 code.  It treats rotation the same in all
    **    cases, because that is correct.  Let the engine rotate the unit
    **    line, and take what it gives you.  Simplicity is beauty.
    **  ptr B723271 -- lotus freelance
    **
    **                                                                   
    */
    if ((pXformData->xform.fxM12) || ((!pXformData->xform.fxM12) &&
       (pXformData->xform.fxM11 == MAKEFIXED(-1, 0))))
    {
      POINTL ptlRot[2];

      pDDC->DCState.LogPosition = PtlTemp[0];

      /*
      ** create a "unit line" to use for the rotation transformation
      ** test.  the value of one-hundred (100) was chosen as it
      ** seemed to provide sufficient accuracy, without risking
      ** coordinate overflow.
      */
      ptlRot[0].x = 0;
      ptlRot[0].y = 0;
      ptlRot[1].x = 100;
      ptlRot[1].y = 0;

      /*
      ** tranform the unit line using the current XformData
      */
      if ((*daConvertWithMatrix)(hDC, (PPOINTL)ptlRot, 2L,
                                 (PXFORM)&pXformData->xform, pDDC,
                                 NGreConvertWithMatrix))
      {
        /*
        **      Can set the rotation - the difference between the
        **  start and finish of the "line" in ptlRot is due to
        **  the rotation component of the transform.
        **
        **
        **  The run - x component
        */
        pDDC->DCState.ptlRotation.x = ptlRot[1].x-ptlRot[0].x;

        /*
        **  The rise - y component
        */
        pDDC->DCState.ptlRotation.y = ptlRot[1].y-ptlRot[0].y;

        /*
        **  Indicate rotation is present.
        */
        pDDC->DCState.usXfData |= XF_ROTATION;
      }
    }
    else
    {
      /*
      **  Indicate rotation isn't present...
      */
      pDDC->DCState.usXfData &= ~XF_ROTATION;
    }

    /*
    **  Save the current matrix so we can avoid repeating the
    **  above process when called again with the same matrix.
    */
    pDDC->DCState.xfLastMatrix = pXformData->xform;

    /*
    ** Viewer shear is not currently addressed, as it is
    ** considered an abberant behavior on the users part,
    ** highly unlikely, and we haven't enough information
    ** to accurately isolate it.
    **
    ** Check if this is a scale only transform.  If so,  remember it,
    ** since it means that some of the plotter's hardware commands can
    ** be used,  rather than calling the engine or simulating ourself.
    */
    pDDC->DCState.usXfData |= XF_SCALE_ONLY;

    if (pXformData->xform.fxM12 || pXformData->xform.fxM21)
    {
      /*
      **  Some other aspects of transform besides scaling!
      */
      pDDC->DCState.usXfData &= ~XF_SCALE_ONLY;
    }
  }

  LeaveDriver(pDDC);
  return  lResult;
}

/***************************************************************************
 *
 * FUNCTION NAME = get_angle_xform
 *
 * DESCRIPTION   = Generate a rotation matrix from character direction
 *                 data.
 *
 * INPUT         = pptlAngle - input angle
 *
 * OUTPUT        = pxfmRot   - rotational matrix output
 *
 * RETURN-NORMAL = TRUE if rotation is applied,  results in pxfmRot.
 * RETURN-ERROR  = FALSE if there is no rotation.
 *         
 *                                                                    
 *                                                                     
 *
 **************************************************************************/

BOOL  get_angle_xform ( PXFORM pxfmRot,/* Rotational matrix output goes here                         */
                        PPOINTL pptlAngle)/* Input angle is received here                     */
/*
** pxfmRot    Rotational matrix output goes here
** pptlAngle  Input angle is received here
*/
{
  /*
  **  Examine the character rotation information,  and generate a rotation
  **  transform if so required.  Returns TRUE if there is a rotation,
  **  else FALSE indicating no rotation (NOT related to character direction).
  **    Output Matrix:
  **          M11 = anglex / sqrt( yangle ** 2 + xangle ** 2 )
  **          M12 = angley / sqrt( yangle ** 2 + xangle ** 2 )
  **          M21 = -angley / sqrt( yangle ** 2 + xangle ** 2 )
  **          M22 = anglex / sqrt( yangle ** 2 + xangle ** 2 )
  */
  LONG lSqrt;
  LONG lX,lY;                          /* Local copies of the angle data    */
  BOOL bXSign,bYSign;                  /* Remember the signs of the values  */

  lX = pptlAngle->x;
  lY = pptlAngle->y;

  while (lX < -65536 || lX > 65535 || lY < -65536 || lY > 65535)
  {
    lX = lX/2;                  /* Prevent overflow of the square function  */
    lY = lY/2;
  }

  /*
  **    Verify that the angles are reasonable - i.e. both are not zero.
  */
  if (lX == 0)
  {

    if (lY == 0)
    {
      GplErrSetError(PMERR_INV_CHAR_ANGLE_ATTR);
      return  FALSE;
    }
  }

  /*
  **   Break out cases where we can avoid loss in precision
  **   Corrected the 90 and 270 cases to assume CCW rotation
  **                                                                     
  */
  if ((!lX) && (lY))
  {                                        /* 90 or 270                     */
    pxfmRot->fxM11 = 0L;
    pxfmRot->fxM22 = 0L;

    if (lY < 0)
    {                                      /* must be 270                   */
      pxfmRot->fxM12 = MAKEFIXED(-1, 0);   /* (sin)                         */
      pxfmRot->fxM21 = MAKEFIXED(1, 0);    /* -(sin)                        */
    }
    else
    {                                      /* 90                            */
      pxfmRot->fxM12 = MAKEFIXED(1, 0);    /* (sin)                         */
      pxfmRot->fxM21 = MAKEFIXED(-1, 0);   /* -(sin)                        */
    }
    pxfmRot->lM41 = 0L;
    pxfmRot->lM42 = 0L;
    return (TRUE);
  }

  if ((!lY) && (lX))
  {                                        /* 0 or 180                      */
    pxfmRot->fxM12 = 0L;
    pxfmRot->fxM21 = 0L;

    if (lX < 0)
    {                                      /* 180                           */
      pxfmRot->fxM11 = pxfmRot->fxM22 = MAKEFIXED(-1, 0);
    }
    else
    {                                      /* 0                             */
      pxfmRot->fxM11 = pxfmRot->fxM22 = MAKEFIXED(1, 0);
    }
    pxfmRot->lM41 = 0L;
    pxfmRot->lM42 = 0L;
    return (TRUE);
  }

  /*
  **  We are now in business - onto the real work.
  */
  lSqrt = f_lsqrt((lX *lX)+(lY *lY));

  if (lX < 0)
  {
    lX = -lX;
    bXSign = TRUE;
  }
  else
    bXSign = FALSE;

  if (lY < 0)
  {
    lY = -lY;
    bYSign = TRUE;
  }
  else
    bYSign = FALSE;
  pxfmRot->fxM11 = f_ldivf(lX, lSqrt);
  pxfmRot->fxM12 = f_ldivf(lY, lSqrt);
  pxfmRot->fxM21 = f_ldivf(lY, lSqrt);
  pxfmRot->fxM22 = f_ldivf(lX, lSqrt);
  pxfmRot->lM41 = 0L;
  pxfmRot->lM42 = 0L;                      /* No translation                */

  /*
  ** Adjust for the different quadrants,  as required
  **
  ** the rotation transformation matrix takes the following
  ** form for counter-clockwise rotation (which is how a
  ** character angle is expressed
  **
  **      Ú                   ¿
  **      ³ cos    sin     0  ³
  **      ³ -sin   cos     0  ³
  **      ³ 0      0       1  ³
  **      À                   Ù
  **
  ** the previous implementation was assuming CW rotation, which
  ** is incorrect for a character angle.  fixed for
  ** ptr B723271 -- lotus freelance
  **
  **                                                                     
  */
  if (bXSign)
  {
    /*
    ** if the x value was originally signed, then
    ** the cosine values must both be re-signed.
    */
    pxfmRot->fxM11 = -pxfmRot->fxM11;
    pxfmRot->fxM22 = -pxfmRot->fxM22;
  }

  if (bYSign)
  {
    /*
    ** if the original y value was less than zero, then
    ** the resulting sine values need to be re-negated
    */
    pxfmRot->fxM12 = -pxfmRot->fxM12;

    /*
    ** fxM21 is the -sin() and is therefore already in its
    ** correct state
    */
  }
  else
  {
    /*
    ** if the original y value was positive, ensure that
    ** fxM21 is a negative value
    */
    if (pxfmRot->fxM21 > 0)
    {
      pxfmRot->fxM21 = -(pxfmRot->fxM21); /* negate it                      */
    }
  }
  return  TRUE;
}

/***************************************************************************
 *
 * FUNCTION NAME = f_lsqrt
 *
 * DESCRIPTION   = Returns the square root of the unsigned long
 *                 argument.  The value is returned as a FIXED number!
 *                 Be warned.
 *
 * INPUT         = ulIn - integer to square root
 *
 * OUTPUT        = FIXED square root
 *
 * RETURN-NORMAL = NONE
 * RETURN-ERROR  = NONE
 *         
 *                                                                      
 *
 **************************************************************************/

ULONG  f_lsqrt ( ULONG ulIn )
{

  /*
  **    Steps in the calculation are:
  **  1/ Normalize the input value (set top bit to 1)
  **  2/ Approximate the square root using Newton's method
  **  3/ Remove the effects of the normalization.
  */

  /*
  **  Number of bits shifted for normalisation
  */
  INT cNormal;


  if (!(ulIn))
    return 0L;

  /*
  **     Perform the normalization.  The reason for this is to maintain
  **  maximum precision.
  */
  cNormal = 0;

  while (!(ulIn&0xff000000))
  {
    ulIn <<= 8;                        /* 8 bits at a time                  */
    cNormal += 8;
  }

  /*
  **  Now for the remaining bits
  */
  while (!(ulIn&0x80000000))
  {
    ulIn <<= 1;
    cNormal++;
  }

  if (cNormal&0x1)
  {
    /*
    **   NOTE:  the shift was by an odd number of bits.  This is not
    ** acceptable, since later we need to shift right by HALF the number
    ** of places that we shifted left above.  So,  reduce the count and
    ** shift the value one place to the right.
    */
    --cNormal;
    ulIn >>= 1;
  }

  /*
  **     Value is normalised to maximise precision.  Now proceed to use
  **  Newton's method to generate the square root.
  */
  ulIn = newton_f_lsqrt(ulIn);

  /*
  **    Now have a fixed version of the square root.   Unnormalize it to
  **  compensate for what was done above.
  */
  cNormal = cNormal/2-16;

  if (cNormal > 0)
    ulIn >>= cNormal;
  else
  {
    if (cNormal < 0)
      ulIn <<= -cNormal;
  }
  return  ulIn;
}

/***************************************************************************
 *
 * FUNCTION NAME = newton f_lsqrt
 *
 * DESCRIPTION   = Evaluate a square root using Newton's approximation
 *                 method.  Uses a fixed number of guesses.
 *
 * INPUT         = ulIn - integer to square root
 *
 * OUTPUT        = A fixed point representation of the square root.
 *
 * RETURN-NORMAL = NONE
 * RETURN-ERROR  = NONE
 *         
 *
 **************************************************************************/

ULONG  newton_f_lsqrt(ULONG ulIn)
{

  /*
  **    Returns an approximation to the square root of ulIn.  Value is
  **  as a fixed point number (16 bits of integer, 16 bits of fraction).
  */
  ULONG ulGuess;


  if (ulIn == 0)
    return 0;                          /* Easy one!                         */
  ulGuess = (ulIn >> 16)+0x4000;

  /*
  ** Next guess = (In / guess + guess) / 2
  */
  ulGuess = (ulIn/ulGuess+ulGuess)/2;

  /*
  ** And one more time for good measure
  */
  ulGuess = (ulIn/ulGuess+ulGuess)/2;

  /*
  ** And one more time for good measure
  */
  ulGuess = (ulIn/ulGuess+ulGuess)/2;
  return  ulGuess;
}

/***************************************************************************
 *
 * FUNCTION NAME = f_ldivf
 *
 * DESCRIPTION   = Divide long by fixed,  returning a fixed.
 *
 * LIMITATIONS
 *                 This function was written with a specific application
 *                 in mind, namely generating the sine/cosine values in
 *                 a triangle.  As such, it may have restrictions not
 *                 suitable for general use.  The result is assumed to be
 *                 less than or equal to 1.0000
 *
 * INPUT         = ulIn  - integer numerator
 *                 uldiv - FIXED divisor
 *
 * OUTPUT        =
 *
 * RETURN-NORMAL = NONE
 * RETURN-ERROR  = NONE
 *         
 *
 **************************************************************************/

ULONG  f_ldivf ( ULONG ulIn,  /* Integer numerator               */
                 ULONG uldiv) /* FIXED divisor                   */
{

  /*
  **    Returns ulIn divided by ulDiv.  Assumes all values are positive.
  */
  INT cCount;                          /* Record number of shifts involved  */


  if (ulIn)
  {

    /*
    **  Non zero,  so off to work.   The test for zero is required
    **  because ulIn is shifted left until the MSB is 1 - this maximises
    **  precision in the result,  but will fail if the number is zero!
    */
    cCount = 0;

    /*
    **      Convert the denominator to a fixed point value.  Note
    **  that it is assumed this will not cause overflow,  as these
    **  values are tested by our caller - specific to angle_xform().
    */
    ulIn <<= 16;

    /*
    **  Shift to get the MSB a 1!
    */
    while (!(ulIn&0x80000000))
    {
      ulIn <<= 1;
      ++cCount;
    }

    /*
    **      We now want to divide a shifted fixed by a fixed,  and still
    **  produce a fixed.   To do this,  we turn the divisor into a long
    **  by shifting it right 16 places.  HOWEVER,  since we have shifted
    **  the numerator LEFT by cCount bits,  we need to shift the
    **  denominator by (16 - cCount) bits to compensate for the shift
    **  of the numerator.  This way,  we maintain maximum precision.
    */

    return  ulIn/(uldiv >> (16-cCount));
  }

  else

    return 0L;                         /* Speedy                            */
}

