TxDuino: Configuration Builder in WxWidgets


In relation to the TxDuino (somewhat documented here), I wrote a C++ class to ease the interfacing with the device. That code was posted here. In this post, I am sharing a program that can be used to easily discover and save the configuration of the saturation points of the different actuators on a vehicle. The program is written using the wxWidgets toolkit so should be compatible with any wx- supported operating system. Of course, at this point in time I still haven't updated the TxDuino class with linux implementations yet, but hopefully that will happen soon.

The saturation points are the raw commands that correspond to the extreme edges of the actuators range. For instance, if the actuator is a servo, then the minimum saturation point is the minimum command that the servo can actually achieve. If a command lower than that is sent, then the servo will constantly jitter as it cannot actually achieve that rotation. By setting the saturation points, normalized commands (percentages) can be sent to the device using the TxDuino::setPercent() function.

A screenshot of the program in action is below:

Actuator Configuration Tester

You can connect to the device by entering the device name in the text field at the top and clicking "Connect" (or pressing enter). Once the device is connected (success or failure will be reported in the console at the bottom), then commands can start being sent. You can start sending commands by clicking "Continue" or stop at any time by clicking "Pause". Note that "pausing" won't stop your vehicle from doing something retarded, the vehicle's actuators will just be stuck in the last commanded state. If the program is not paused, then every time you change one of the values in the spin controls, then that new value is sent as the current command.

To discover the value of the minimum saturation point of one of the actuators, start by clicking on the spin control for that channel, setting it to zero, and hitting enter (this makes the control think it has changed whether it has or not so the value is sent to that channel). If the actuator is saturated (on a servo, if it's making noise and/or jittering but not correcting itself) then raise the value by one unit. Keep raising the value until the actuator is no longer saturated. Repeat for the max saturation point and the neutral point. The neutral point is the point corresponding to 0%. For a servo this is most likely near the center of it's range. For a speed controller, it should be the same as the minimum saturation.

A windows binary of this program can be downloaded here: Tester .
The source can be downloaded here: TxDuino Actuator Configuration Tester.

TxDuino: C++ classes (windows)


In order to make interfacing with the TxDuino as easy as possible, I wrote a little c++ class to take care of all the heavy lifting. You connect to a TxDuino device by creating a TxDuino object. The constructor accepts a device name. In windows that looks like "\.COM8" (my arduino is installed on COM port 8). In linux it'll look something like "/dev/ttyS2" except that I haven't gotten around to implementing the linux part of this class yet.

The class, along with the supporting classes and test program are included here: TxDuino Class Source.

/**
 *  file       TxDuino.h
 *  date:      Oct 27, 2009
 *  brief:
 *
 *  detail:
 */

#ifndef CTXDUINO_H_
#define CTXDUINO_H_


#include 
#include 

#include "types.h"


namespace txduino
{


/**
 *  brief  operating system dependendant implementation structure
 */
struct STxDuinoImpl;

/**
 *  brief Class encapsulating the interface / API for communicating with one
 *  TxDuino device. The TxDuino sends a standard RC PPM signal encoding
 *  commands for up to 8 PWM channels (servos, engine controllers).
 *
 *  The magnitude of each channel is divided into 250 discrete segments.
 *  Exactly how those segments are interpreted by the actuators is determined
 *  by the configuration of this object.
 *
 *  note: channels are zero indexed
 *
 *  todo   add a constructor that accepts a csv file with the saturations and
 *          neutral points that can be generated from the test program
 */
class CTxDuino
{
    private:
        STxDuinoImpl*   m_osInfo;   /// pointer to os-dependant implementation
        std::string     m_devName;  /// system device name,
                                    /// i.e. "\.COM2″, "/dev/tty2″

        u8          m_chan      [9];    /// value for each channel [0,250]
        u8          m_neutral   [8];    /// the "center" for each actuator
        u8          m_minsat    [8];    /// the minimum value for each channel
        u8          m_maxsat    [8];    /// the maximum value for each channel


    public:
        /**
         *  brief  Constructs a new TxDuino object serving as an interface
         *          into on particular TxDuino device.
         *  param  strDevice   device string, i.e. "\.COM2″, "/dev/tty2″
         */
        CTxDuino( std::string strDevice );


        /**
         *  brief  cleans up OS resources reserved for this serial connection
         */
        virtual ~CTxDuino();


        /**
         *  brief  sends the current channel definitions to the device
         */
        void send();


        /**
         *  brief  sets the value of an actuator as a percent of it's viable
         *          range.
         *  param  chan        the channel to set
         *  param  percent     -100% < percent < 100%; value to set channel to
         */
        void setPercent( s32 chan, f64 percent );


        /**
         *  brief  returns the percent value the indicated actuator is set to,
         *          note: if the neutral point is equal to one of the saturation
         *          points this value may be unreliable
         *  param  chan        the channel to get
         *  return percent value of actuator on indicated channel
         */
        f64 getPercent( s32 chan );


        /**
         *  brief  sets the raw value of the actuator pulse width
         *  param  chan        the channel to set
         *  param  value      0 < value < 250; value to set channel to
         */
        void setRaw( s32 chan, u8 value );


        /**
         *  brief  returns the raw value the indicated actuator is set to
         *  param  chan        the channel to get
         *  return a value between 0 and 250 indicated the pulse length for
         *          that actuator (multiply by 4us and add 700us to get the
         *          actual pulse length)
         */
        u8 getRaw( s32 chan );




        /**
         *  brief  sets the raw value of the actuator pulse width corresponding
         *          to a neutral state of that actuator
         *  param  chan        the channel to set
         *  param  value      0 < value < 250; value to set channel to
         */
        void setNeutral( s32 chan, u8 value );


        /**
         *  brief  returns the raw value corresponding to a neutral state of
         *          the indicated actuator
         *  param  chan        the channel to get
         *  return a value between 0 and 250 indicated the pulse length for
         *          that actuator (multiply by 4us and add 700us to get the
         *          actual pulse length)
         */
        u8 getNeutral( s32 chan );




        /**
         *  brief  sets the raw value of the actuator pulse width corresponding
         *          to the minimum state of the indicated actuator
         *  param  chan        the channel to set
         *  param  value      0 < value < 250; value to set channel to
         */
        void setMinSat( s32 chan, u8 value );


        /**
         *  brief  returns the raw value corresponding to a minimum state of
         *          the indicated actuator
         *  param  chan        the channel to get
         *  return a value between 0 and 250 indicated the pulse length for
         *          that actuator (multiply by 4us and add 700us to get the
         *          actual pulse length)
         */
        u8 getMinSat( s32 chan );




        /**
         *  brief  sets the raw value of the actuator pulse width corresponding
         *          to the maximum state of the indicated actuator
         *  param  chan        the channel to set
         *  param  value      0 < value < 250; value to set channel to
         */
        void setMaxSat( s32 chan, u8 value );


        /**
         *  brief  returns the raw value corresponding to a maximum state of
         *          the indicated actuator
         *  param  chan        the channel to get
         *  return a value between 0 and 250 indicated the pulse length for
         *          that actuator (multiply by 4us and add 700us to get the
         *          actual pulse length)
         */
        u8 getMaxSat( s32 chan );


        /**
         *  brief  return the device name that was used to connect to this
         *          txduino
         */
        std::string getName();

};

}

#endif /* CTXDUINO_H_ */



/**
 *  file       CTxDuino.cpp
 *  date:      Oct 27, 2009
 *  brief:
 *
 *  detail:
 */

#include "CTxDuino.h"
#include "compile.h"

#include 
#include 
#include 
#include 

#include "IllegalArgumentException.h"
#include "IOException.h"

#ifdef TXD_MINGW
#include 
#endif

namespace txduino
{



#ifdef TXD_MINGW
struct STxDuinoImpl
{
    HANDLE hComPort;    /// handle to the opened COM device
};
#endif


#ifdef TXD_LINUX
struct STxDuinoImpl
{
    FILE hSerialFile;
};
#endif






/**
 *  The initial state of the individual actuators is initialized as follows
 *
 *  verbatim
 *      minimum saturation: 0
 *                 neutral: 125
 *      maximum saturation: 250
 *                   value: 125
 *  endverbatim
 *
 *  Note that this is probably not appropriate for your system. Many of the
 *  servos that we've used have minimum saturations around 10 and maximum
 *  saturations around 240. Use the actuator test program to determine what
 *  these values should be.
 */
CTxDuino::CTxDuino( std::string strDevice )
{
    using std::cout;
    using std::endl;
    using std::stringstream;

    // intiialize all the arrays
    for( int i=0; i < 8; i++ )
    {
        m_chan      [i] = 125;
        m_minsat    [i] = 0;
        m_maxsat    [i] = 250;
        m_neutral   [i] = 125;
    }

    // stop byte
    m_chan[8]   = 0xFF;

    // initialize the OS dependant information
    m_osInfo    = new STxDuinoImpl();


/* ----------------------------------------------------------------------------
 * Windows Specific Implementation:
 * ---------------------------------------------------------------------------*/

#ifdef TXD_MINGW
    // open the file using the windows API
    m_osInfo->hComPort  =
    CreateFile( strDevice.c_str(),          // file name
        GENERIC_READ | GENERIC_WRITE,       // access mode: read and write
        FILE_SHARE_READ|FILE_SHARE_WRITE,   // (sharing)
        NULL,                               // (security) 0: none
        OPEN_EXISTING,                      // (creation) i.e. don't make it
        0,                                  // (overlapped operation)
        NULL);                              // no template file

    // check to make sure the file open succeeded
    if( m_osInfo->hComPort == INVALID_HANDLE_VALUE )
    {
        stringstream message( stringstream::in | stringstream::out );
        message << "Invalid Device Name: " << strDevice;
        throw IllegalArgumentException(message.str());
    }

    // get the current settings on the com port
    DCB dcb;
    GetCommState( m_osInfo->hComPort, &dcb );

    // change the settings, the TxDuino uses a BAUD rate of 9600
    dcb.fBinary     =   1;
    dcb.BaudRate    =   CBR_9600;
    dcb.Parity      =   NOPARITY;
    dcb.ByteSize    =   8;
    dcb.StopBits    =   ONESTOPBIT;

    // set the new settings for the port
    SetCommState( m_osInfo->hComPort, &dcb );
#endif


}



CTxDuino::~CTxDuino()
{

/* ---------------------------------------------------
 * Windows Specific Implementation:
 * --------------------------------------------------*/

#ifdef TXD_MINGW
    // close the device if it's open
    if( m_osInfo->hComPort != INVALID_HANDLE_VALUE )
        CloseHandle( m_osInfo->hComPort );
#endif

delete m_osInfo;

}



void CTxDuino::send()
{
    using std::stringstream;


    // ensure that the last byte of the packet is the stop byte
    m_chan[8] = 0xFF;


/* ---------------------------------------------------
 * Windows Specific Implementation:
 * --------------------------------------------------*/

#ifdef TXD_MINGW
    DWORD bytesWritten;

    BOOL retVal =
    WriteFile(  m_osInfo->hComPort, // output handle
                m_chan,             // buffer of bytes to send
                9,                  // number of bytes to send from buffer
                &bytesWritten,      // pointer to a word that receives number of
                                    // bytes written
                NULL);              // pointer to an OVERLAPPED struct

    if( bytesWritten != 9 )
    {
        stringstream message( stringstream::in | stringstream::out );
        message << "Bytes written to device less than expected: "
                << bytesWritten << ", expecting 9";
        throw IOException(message.str());
    }

    if( retVal == 0 )
    {
        stringstream message( stringstream::in | stringstream::out );
        message << "Writing to device failed; error code: " << GetLastError();
        throw IOException(message.str());
    }
#endif

}



/**
 *  If the percent is positive, the raw value is calculated as follows
 *
 *  verbatim
 *      raw = neutral + (max - neutral) * percent
 *  endverbatim
 *
 *  if the percent is negative, the raw value is calculated as follows
 *
 *  verbatim
 *      raw = neutral - (neutral - min) * percent
 *  endverbatim
 */
void CTxDuino::setPercent( s32 chan, f64 percent )
{
    using std::stringstream;

    if(chan < 0 || chan > 7)
    {
        stringstream message( stringstream::in | stringstream::out );
        message << "Invalid channel number: " << chan << HERE
                << "valid channels are 0-7";
    }

    if(percent > 0)
    {
        m_chan[chan] = (u8) (m_neutral[chan] +
                    ( m_maxsat[chan] - m_neutral[chan]) * percent );
    }

    else
    {
        m_chan[chan] = (u8) (m_neutral[chan] -
                    ( m_minsat[chan] - m_neutral[chan]) * percent );
    }
}



/**
 *  If the value is strictly less than the neutral value then the percent value
 *  is calculated by
 *
 *  verbatim
 *      percent = -( neutral - value ) / (neutral - min);
 *  endverbatim
 *
 *  If the value is strictly greater than the neutral value then the percent
 *  value is calculated by
 *
 *  verbatim
 *      percent = ( value - neutral ) / (max - neutral);
 *  endverbatim
 *
 *  If the value is equal to the neutral value then the percent value is zero.
 */
f64 CTxDuino::getPercent( s32 chan )
{
    using std::stringstream;

    if(chan < 0 || chan > 7)
    {
        stringstream message( stringstream::in | stringstream::out );
        message << "Invalid channel number: " << chan << HERE
                << "valid channels are 0-7";
    }

    if( m_chan[chan] < m_neutral[chan] )
    {
        return (double)(-(m_neutral[chan] - m_chan[chan])) /
                    (m_neutral[chan] - m_minsat[chan]);
    }

    else if( m_chan[chan] > m_neutral[chan] )
    {
        return (double)(m_chan[chan] - m_neutral[chan]) /
                    (m_maxsat[chan] - m_neutral[chan]);
    }

    else
    {
        return 0.0;
    }
}



void CTxDuino::setRaw( s32 chan, u8 value )
{
    using std::stringstream;

    if(chan < 0 || chan > 7)
    {
        stringstream message( stringstream::in | stringstream::out );
        message << "Invalid channel number: " << chan << HERE
                << "valid channels are 0-7";
    }

    if(value > 250)
    {
        stringstream message( stringstream::in | stringstream::out );
        message << "Invalid value: " << value << HERE
                << "valid channels are 0-250";
    }

    m_chan[chan] = value;
}



u8 CTxDuino::getRaw( s32 chan )
{
    using std::stringstream;

    if(chan < 0 || chan > 7)
    {
        stringstream message( stringstream::in | stringstream::out );
        message << "Invalid channel number: " << chan << HERE
                << "valid channels are 0-7";
    }

    return m_chan[chan];
}



void CTxDuino::setNeutral( s32 chan, u8 value )
{
    using std::stringstream;

    if(chan < 0 || chan > 7)
    {
        stringstream message( stringstream::in | stringstream::out );
        message << "Invalid channel number: " << chan << HERE
                << "valid channels are 0-7";
    }

    if(value > 250)
    {
        stringstream message( stringstream::in | stringstream::out );
        message << "Invalid value: " << value << HERE
                << "valid channels are 0-250";
    }

    m_neutral[chan] = value;
}



u8 CTxDuino::getNeutral( s32 chan )
{
    using std::stringstream;

    if(chan < 0 || chan > 7)
    {
        stringstream message( stringstream::in | stringstream::out );
        message << "Invalid channel number: " << chan << HERE
                << "valid channels are 0-7";
    }

    return m_neutral[chan];
}



void CTxDuino::setMinSat( s32 chan, u8 value )
{
    using std::stringstream;

    if(chan < 0 || chan > 7)
    {
        stringstream message( stringstream::in | stringstream::out );
        message << "Invalid channel number: " << chan << HERE
                << "valid channels are 0-7";
    }

    if(value > 250)
    {
        stringstream message( stringstream::in | stringstream::out );
        message << "Invalid value: " << value << HERE
                << "valid channels are 0-250";
    }

    m_minsat[chan] = value;
}



u8 CTxDuino::getMinSat( s32 chan )
{
    using std::stringstream;

    if(chan < 0 || chan > 7)
    {
        stringstream message( stringstream::in | stringstream::out );
        message << "Invalid channel number: " << chan << HERE
                << "valid channels are 0-7";
    }

    return m_minsat[chan];
}



void CTxDuino::setMaxSat( s32 chan, u8 value )
{
    using std::stringstream;

    if(chan < 0 || chan > 7)
    {
        stringstream message( stringstream::in | stringstream::out );
        message << "Invalid channel number: " << chan << HERE
                << "valid channels are 0-7";
    }

    if(value > 250)
    {
        stringstream message( stringstream::in | stringstream::out );
        message << "Invalid value: " << value << HERE
                << "valid channels are 0-250";
    }

    m_maxsat[chan] = value;
}



u8 CTxDuino::getMaxSat( s32 chan )
{
    using std::stringstream;

    if(chan < 0 || chan > 7)
    {
        stringstream message( stringstream::in | stringstream::out );
        message << "Invalid channel number: " << chan << HERE
                << "valid channels are 0-7";
    }

    return m_maxsat[chan];
}



std::string CTxDuino::getName()
{
    return m_devName;
}

}

A re-write of the serialTest.exe program that sends sinusoidal commands to the plane using this new class demonstrates its use.

/**
 *  file        serialTest.cpp
 *  date:      Oct 27, 2009
 *  brief:
 *
 *  detail:
 *  This is a simple test program that demonstrates how to connect to and
 *  write commands to the arduino transmitter interface using windows.
 *
 *  the TxDuino is an interface into the futaba FP-TP-FM transmitter module,
 *  which accepts an RC PPM input. This signal contains a maximum of 8
 *  servo channels.
 */

#include 
#include 
#include 

#include "CTxDuino.h"
#include "constants.h"

using namespace std;
using namespace txduino;

int main( int argc, char** argv )
{
    // check to ensure that the command line included a device to open
    if( argc < 2 )
    {
        cout << "Usage: serialTest.exe [Device Name]n"
                "   where [Device Name] is the name of the COM port file onn "
                "   windows (i.e. \\.\COM8), or the name of the serialn "
                "   device on *nix (i.e. /dev/tty8)n" << endl;

        return -1;
    }

    // grab a pointer to the device to open
    char*       strDevName = argv[1];

    // create the txduino device
    CTxDuino tx(strDevName);

    // send a sinusoidal input on all channels (except for channel 3, which is
    // usually the throttle) for 10 seconds
    for(int i=0; i < 1000; i++)
    {
        for(int j=0; j < 8; j++)
            tx.setRaw(j, (unsigned char)
                            (125.0 + 75.0 * sin( 2.0 * PI * i / 100.0 )) );

        tx.setRaw(2, 0);

        for(int j=0; j < 8; j++)
            cout << setw(3) << (int)tx.getRaw(j) << " | ";
        cout << endl;

        tx.send();

#ifdef TXD_MIGNW
        Sleep(1);
#endif

#ifdef TXD_LINUX
        usleep(0.001);
#endif
    }


    return 0;
}

TxDuino: Custom PC Controlled RC Transmitter


In our lab we have some very light weight foam aircraft that we do control work with. In order to keep the weight down, we use micro RC receivers that generate PWM signals for each of the servos that deflect one of three control surfaces: ailerons, elevator, and rudder. Note that both ailerons are linked and are together controlled by a single servo (a common setup for these lightweight foam planes). A PWM signal is also sent to control the throttle of the propeller. The glowing little spheres are reflective features to aid in motion capture which we use for estimation of the vehicle's state.

One of our 'lil baby airplanes

The trick, however, is to get our control codes, running on a PC somewhere in the room, to actually change the deflection of the airplane's control surfaces, or the throttle of the propeller.

The TxDuino project aims to create a single device that allows a PC to control an RC vehicle. At the time of starting this project, the only option for doing this was a device from here that could be plugged into the trainer port of a transmitter. An arduino has been used to achieve this same setup as described here. Going through the handheld transmitter, however, is a huge hassle for a number of reasons so I set out to create a simple single device that could accomplish this. As it turns out, so did the guys at Tom's RC, their solution being this guy here. Anyway, since the work is done, I'll share the results.

The easiest way that I figured to pull this off is to work with a modular RC transmitter, and interface directly with the module. For instance this is a transmitter

Futaba Modular Transmitter

And this is the module that plugs into the back of the transmitter

the FP-TP-FM Transmitter Module

The handheld part creates some kind of signal that it sends to the module, which does something, and then sends something else to the antenna. If I can figure out how to recreate that signal, then the "somethings" aren't very important. The modules are cheap enough and easy enough to find that there isn't much of a need to re-engineer that part.

To begin with, I needed to do some investigation. There are 5 pins on the back of the transmitter. It was easy enough by tracing the module circuit to find the source, ground, and RF out (antenna) pins, but what about the other two? I was hoping one of them was the control signal in a not too esoteric form, and I had no idea what the last would be. I managed to find the control signal with a little probing. Here is a little video showing the transmitter (Tx) hooked up to the module, via some wires with alligator clips which made it easy to investigate things.

After figuring our a little of what was going on, I found this forum post that had a pretty good discussion of the pinouts on this particular module. In particular, this post I found to be accurate and helpful.

Basically, looking at the back of the transmitter (or the top of the module, when the cover is off) the pins are

    1   2   3   4   5
    .   .   .   .   .

 left side
 #1 PPM
 #2 V+
 #3 RF check
 #4 ground
 #5 RF out
 right side

Of particular importance is the note on that post that the hand-held part is actually pulling the line low to make the pulses, and that the PPM line is at 2.5 volts.

I found a description of the PPM signal here. While the measurements I made using the oscilloscope did not match the information on that page exactly (I measured a 4us pulse between channels), using the signal description from that page has had good results. The details of the PPM frame are repeated here.

  • complete PPM-Frame has a length of 22.5 ms.
  • consists of an overlong start segment combined with 8 analog channels.
  • stop pulse of 0.3 ms follows every channel's impulse.
  • The length of a channel's impulse ranges from 0.7ms to 1.7 ms

Below is a video of the PPM signal from the transmitter. You can see the impulse length for each channel changes as I manipulate the control sticks on the handheld part of the transmitter. I apologize for the poor light sensing on my crappy little digital camera. Also, note that channel 5 jumps up and down because it is a flip switch (see the photo above) and only has two possible positions.

Ok, now that I know what is happening, I need a way to recreate that PPM signal using some kind of device that I can connect to my PC. Hey, the arduino is easy to develop with, comes with a USB controller, and has a nice software library that handles the PC interface stack. It's also built around an Atmel microcontroller (my decimila specifically uses a ATmega168L) which has some pretty sophisticated features. It runs at 16Mhz, has three timers that can be used for generating interrupts, and comes with a (relatively) standard C / C++ compiler for writing code. Furthermore, the arduino's bootloader means that I don't even have to worry about messing with a programmer to get something up and running.

I wont go into all the gory details of the TxDuino design, but I will point out some of the relevant considerations.

The 168L has three timers. Timer1 is a 16-bit counter, Timer0 and Timer2 are 8-bit counters. Timers 0 and 1 share a prescaler, Timer2 has it's own. The prescaler is a component that reduces the update rate of the timer, so we can scale the 16Mhz system clock to some fraction of it's rate. Using an 8-bit counter would be extremely limiting because scaling the clock down enough that the counter capacity is sufficient would provide a counter resolution that is too coarse for the fine control we need to generate our waveform. This means I want to use Timer1. However, Timer0 is used by the serial stack. I want to be able to send commands from the PC to the TxDuino which will require use of the serial interface. Thus, the serial baud rate will dictate the prescale factor used by Timer0, and thus used by Timer1… meaning that my choices are limited to whatever the serial stack allows.

A simple arduino test program, however, quells this dilemma.

 int inByte = 0;         // incoming serial byte

 //define the Timer1 overflow interrupt service routine
 ISR(TIMER1_OVF_vect) 
 {
    // print out the current value of the timer control register B which
    // contains the timer prescale value, set when the serial BAUD rate 
    // is defined
    Serial.print("Timer1: ");
    Serial.print(TCCR1B, HEX);
    Serial.print("rn");
 }


 void setup()
 {
   // set Timer1 to operate in "normal" mode
   TCCR1A = 0;

   // activate Timer1 overflow interrupt 
   TIMSK1 = 1 << TOIE1;

   // start serial port at 9600 bps:
   Serial.begin(9600);
 }

 void loop()
 {
   // if we get a valid byte, read analog ins:
   if (Serial.available() > 0) {
     // get incoming byte:
     inByte = Serial.read();
     Serial.print(inByte, BYTE);
     Serial.print("rn");
   }
 }

This program sets the serial baud rate to 9600 bits/sec, and then prints out the value of the Timer1 control register. A subset of the bits in that register indicate what the prescaler is set to, and thus what the prescale factor is. Running this program, and looking at the datasheet, I was able to determine that the prescale factor is set to 64. This means that the timer resolution for Timer1 is

$$ \frac{ 64 \frac{ \text{system ticks}}{ \text{counter tick}} } {16,000,000 \frac{ \text{system ticks}}{ \text{second}}} = 0.000004 s = 4 \mu s$$

And the 16-bit timer has a range of zero to

$$ (2^{16} - 1) \cdot (4 \mu s) = 65535 \cdot 4 \mu s = 262140 \mu s = 262+ ms $$

Since our frame length is 22.5 ms, we have plenty of range to work with. Also, the range on the individual servo channels is 1000 us, so that gives us 1000/4 = 250 possible steps for each channel (251 possible discrete positions). That gives us a resolution of better than +/- 0.5% for each channel, which is quite acceptable.

Also, since the PPM signal pulls the voltage on the line low, we cannot just use digitalWrite() like we usually do when writing to a pin with the arduino. Instead, we need a way to drive the voltage to zero when we want to pull the line low, and a way to basically "do nothing" when we want the line high. We can do this by setting the pin voltage to LOW, and then switching between INPUT and OUTPUT using pinMode()

The last consideration is the serial packet format. I decided to keep things simple. Since each channel can only have 251 possible values, I used an unsigned 8-bit integer ($latex 2^8 - 1 = 255$ possible values for an 8-bit number) for each channel, with values 0-250, and reserving 0xFF (255 in decimal) as a end-of-packet identifier.

Here is a little test code to demonstrate the use of pinMode() to set the voltage of the line high / lo. It listens for either an 'a' or 'b' character over serial, and sets the line high or low depending on what is sent. If you want to play with this, you can use the terminal from within the Arduino software to send the commands.

int inByte      = 0;    // incoming serial byte
int ppmPin      = 2;    // pin to make the waveform on


// to setup the device, we will initialize all the channels, and initialize the
// timer settings
void setup()
{
    // start serial port at 9600 bps, note that this also sets the Timer1
    // prescaler, so we need to ensure that we don't change it later; a test 
    // using a previous program showed that 9600 baud set the prescaler to 64
    // or the Timer1 control register B to 0x03.
    Serial.begin(9600);

    // we configure the signal pin as an input; not that the signal actually is
    // an input but we modify the signal by pulling it low when we need to
    pinMode(ppmPin, INPUT);


    // then we set the inital state of the pin to be low
    digitalWrite(ppmPin, LOW);
}


// the main routine monitors the serial port, and on receipt of a completed 
// packet will recalculate the compare points for the specified signal
void loop()
{
    // if there's a byte available on the serial port, then read it in
    if (Serial.available() > 0)
    {
        // get incoming byte
        inByte = Serial.read(); 

        switch(inByte)
        {
            case 'a':
                pinMode(ppmPin, OUTPUT);
                Serial.print("setting pin lown");
                break;

            case 'b':
                pinMode(ppmPin, INPUT);
                Serial.print("setting pin highn");
                break;

            default:
                break;
        }
    }
}

Here is a photo of the arduino hooked up with a resistor from the PPM pin to ground which allows me to look at the signal as I'm testing. The wire from the right side of the resistor to the ground terminal on the arduino is hidden behind the black probe. The other things connected to the breadboard are the transmitter module (the thing in the metal shield), and the back circuit board from the transmitter (the green thing; the pins from that board are plugged into the breadboard). The red thing is the battery for the transmitter.

Arduino Test Setup

Well that's pretty much it as far as background and design. Here is the final code that I used to generate the desired waveform. I wont describe it much here since it's heavily commented already. Basically it uses two compare match interrupts with Timer1. The two interrupts leap-frog each other and signal either a rising edge or a falling edge. The counter is not cleared when the interrupts are serviced, until the last falling edge is serviced.

int inByte      = 0;    // incoming serial byte
int ppmPin      = 2;    // pin to make the waveform on
int chan[8];            // 8 channels encoded as 250 * (percent range)

int compA[9];           // 9 compare points for Timer1
int compB[9];           // 9 compare points for Timer2
int iCompA      = 0;    // compare point index for currently active
int iCompB      = 0;    // compare point index for currently active


// the 0xFF byte will serve as a packet delimeter, we don't want to read
// garbage from the serial so we'll start as soon as we get one of those but
// not before
int bDelimReceived  = 0;

// we don't want to write to any of the compare registers unless we have to 
// so we'll only do so at the end of a frame if that frame involved a change
// in the current signal
int bSignalChanged  = 0;

// index for which channel is next to read in
int iChan = 0;


// define the Timer1 compare match interrupt service routine for the rising 
// edge
ISR(TIMER1_COMPA_vect) 
{
    // to set the line high, we switch the pin mode to "INPUT"
    pinMode(ppmPin, INPUT);


    // if this is not the last rising edge
    if( iCompA < 8 )
    {
        // then we change the compare register to issue an interrupt at the next
        // rising edge, by incrementing the compare point counter
        iCompA++;

        // and then setting the compare register
        OCR1A = compA[iCompA];
    }

    // if this is the last rising edge then we need to set the compare
    // register to be at the time of the rising edge following the start 
    // frame
    else
    {
        // first, if the signal has changed since while writing the last frame
        // then we need to update the compare points
        if( bSignalChanged )
        {
            compA[0] = 3550;
            for(int i = 0; i < 8; i++)
                compA[0] -= chan[i];

            // the end of the first pulse is 75 ticks (300us) after that
            compB[0] = compA[0] + 75;

            // and then the rest of the compare points depend on the value of
            // each channel
            for(int i = 1; i < 9; i++ )
            {
                compA[i] = compB[i-1] + chan[i-1] + 175;
                compB[i] = compA[i] + 75;
            }

            bSignalChanged = 0;
        }

        // then we set the compare point to match the end of the start segment
        iCompA  = 0;
        OCR1A   = compA[0];
    }

}

// define the Timer1 compare match interrupt service routine for the falling 
// edge
ISR(TIMER1_COMPB_vect) 
{
    // to set the pin low we open the collector
    pinMode(ppmPin, OUTPUT);


    // if this is not the last falling edge
    if( iCompB < 8 )
    {
        // then we change the compare register to issue an interrupt at the next
        // rising edge, by incrementing the compare point counter
        iCompB++;

        // and then setting the compare register
        OCR1B = compB[iCompB];
    }

    // if this is the last falling edge
    else
    {
        // then we wrap around to the beginning
        iCompB=0;

        // and set the compare register
        OCR1B = compB[iCompB];

        // and we also need to reset the timer because this is the end of the 
        // frame
        TCNT1 = 0x00;
    }

}






// to setup the device, we will initialize all the channels, and initialize the
// timer settings
void setup()
{
    // initialize all eight channels to be at 125, or neutral
    for(int i=0; i < 8; i++)
        chan[i] = 125;

    // initialize the throttle to be "full left" which is actually "full stop"
    chan[2] = 0;

    // start serial port at 9600 bps, note that this also sets the Timer1
    // prescaler, so we need to ensure that we don't change it later; a test 
    // using a previous program showed that 9600 baud set the prescaler to 64
    // or the Timer1 control register B to 0x03.
    Serial.begin(9600);

    // set Timer1 to operate in "normal" Mode, and set the
    // prescaler to 64; the default mode is one of the PWM modes so we need to
    // set this
    TCCR1A = 0;

    // note: Timer0 and Timer1 share a prescaler, and Timer0 is used by the 
    // serial lines. Therefore, we cannot pick the value of the prescaler, 
    // rather it is determined by the baud rate we want. A simple program has 
    // shown that setting the baud rate to 9600 involves setting the prescaler 
    // to 64

    Serial.print("Timer Control Registers: n");
    Serial.print("   A: 0x");
    Serial.print(TCCR1A, HEX);
    Serial.print("n");
    Serial.print("   B: 0x");
    Serial.print(TCCR1B, HEX);
    Serial.print("n");

    // set the Timer1 output compare register A to 3550 ticks, which corresponds 
    // to 14,200us with the prescale set to 64; the math goes like this: 
    //    with the prescaler set to 64, the Timer1 period is 64/16e6 = 4us
    //    want an interrupt issued after 22.5ms - 8*700us - 9*300us = 14,200us
    //    means we want an interrupt issued after 14,200/4 = 3550 ticks
    compA[0] = 3550;
    OCR1A = compA[0];

    // set Timer1 output compare register B to something 300us/4us = 75 ticks 
    // after the first compare point for timer A
    compB[0] = compA[0] + 75;
    OCR1B = compB[0];

    // we can go ahead and initialize all the other compare points as well, we
    // know that all the channels are set to zero, but I'll do this the long 
    // way because it'll match code that we use elsewhere in this program
    for(int i=1; i < 9; i++)
    {
        compA[i] = compB[i-1] + chan[i-1] + 175;
        compB[i] = compA[i] + 75;
    }

    // activate Timer1 output compare interrupt with compare register A and
    // compare register B by setting the Output Compare Interrupt Enable bits
    TIMSK1 = 1 << OCIE1A | 1 << OCIE1B;


    // we configure the signal pin as an input; not that the signal actually is
    // an input but we modify the signal by an open collector so to make the
    // signal high we enable a 20K pull-up resistor, and to make the signal 
    // low we open the collector
    pinMode(ppmPin, INPUT);


    // then we set the inital state of the pin to be low
    digitalWrite(ppmPin, LOW);
}



// the main routine monitors the serial port, and on receipt of a completed 
// packet will recalculate the compare points for the specified signal
void loop()
{
    // if there's a byte available on the serial port, then read it in
    while (Serial.available() > 0)
    {
        // get incoming byte
        inByte = Serial.read(); 

        // if the byte is a packet delimiter, then start recording the new 
        // packet by resetting the channel indexer
        if( inByte == 0xFF )
        {
            bDelimReceived  = 1;
            iChan           = 0;
        }

        // otherwise, if we've recieved at least one packet delimiter and the
        // current channel index is valid, then record the chanel value; 
        else if(bDelimReceived && iChan < 8)
        {
            // we can save some processing time by only updating the stored 
            // values if the incoming command has changed, so we'll condition
            // the calculations on the fact that the byte is different than
            // what is already stored
            if(inByte != chan[iChan])
            {
                // the protocol is defined to send values between 0 and 250 
                // inclusive, for 250 distinct values; the timing increment is 
                // between 0 and 250 so we simply store the value
                chan[iChan] = inByte;

                // since the signal changed, we need to raise the flag
                bSignalChanged = true;
            }

            // then we increment the channel index so that we're ready for 
            // the next byte that's sent
            iChan++;
        }
    }
}

In order to test this, I wrote this little program to send a sinusoidal input to the TxDuino. Note this code is specific to windows.

/**
 *  file       serialTest.cpp
 *  date:      Oct 27, 2009
 *  brief:
 *
 *  detail:
 *  This is a simple test program that demonstrates how to connect to and
 *  write commands to the arduino transmitter interface using windows.
 *
 *  the TxDuino is an interface into the futaba FP-TP-FM transmitter module,
 *  which accepts an RC PPM input. This signal contains a maximum of 8
 *  servo channels.
 */

#include 
#include 
#include 
#include 

#define PI 3.14159265

int main( int argc, char** argv )
{
    using std::cout;
    using std::endl;
    using std::sin;
    using std::setw;

    // check to ensure that the command line included a device to open
    if( argc < 2 )
    {
        cout << "Usage: serialTest.exe [Device Name]n"
                "   where [Device Name] is the name of the COM port file onn "
                "   windows (i.e. \\.\COM8), or the name of the serialn "
                "   device on *nix (i.e. /dev/tty8)n" << endl;

        return -1;
    }

    // grab a pointer to the device to open
    char*       strDevName = argv[1];

    // open the file using the windows API
    HANDLE      hComPort         =
    CreateFile( strDevName,
        GENERIC_READ | GENERIC_WRITE,       // access mode: read and write
        FILE_SHARE_READ|FILE_SHARE_WRITE,   // (sharing)
        NULL,                               // (security) 0: none
        OPEN_EXISTING,                      // (creation) i.e. don't make it
        0,                                  // (overlapped operation)
        NULL);                              // no template file

    // check to make sure the file open succeeded
    if( hComPort == INVALID_HANDLE_VALUE )
    {
        cout << "FATAL, Invalid device name: " << strDevName << "n" << endl;
        return -2;
    }

    // get the current settings on the com port
    DCB dcb;
    GetCommState( hComPort, &dcb );

    // change the settings, the TxDuino uses a BAUD rate of 9600
    dcb.fBinary     =   1;
    dcb.BaudRate    =   CBR_9600;
    dcb.Parity      =   NOPARITY;
    dcb.ByteSize    =   8;
    dcb.StopBits    =   ONESTOPBIT;

    // set the new settings for the port
    SetCommState( hComPort, &dcb );

    // allocate the data buffer for sending over serial
    unsigned char data[9];

    // set the last byte to be the stop byte
    data[8] = 0xFF;

    // send a sinusoidal input on all channels (except for channel 3, which is
    // usually the throttle) for 10 seconds
    for(int i=0; i < 1000; i++)
    {
        for(int j=0; j < 8; j++)
            data[j] = (unsigned char)(125.0 + 75.0 * sin( 2.0 * PI * i / 100.0 ));

        data[2] = 0;
        data[8] = 0xFF;

        for(int j=0; j < 8; j++)
            cout << setw(3) << (int)data[j] << " | ";
        cout << endl;



        DWORD bytesWritten;

        BOOL retVal =
        WriteFile(  hComPort,       // output handle
                    data,           // buffer of bytes to send
                    9,              // number of bytes to send from buffer
                    &bytesWritten,  // pointer to a word that recieves number of
                                    // bytes written
                    NULL);          // pointer to an OVERLAPPED struct

        Sleep(1);
    }


    // close the device
    CloseHandle( hComPort );

    return 0;
}

The result of running this test program with the probes attached to the arduino is below.

Then I disconnected the line from the handheld part of the transmitter to the module, and connected the arduino's output instead. This is shown below.

Arduino Replacing PPM from Handheld

By connecting the oscilloscope probes like this

Arduino Replacing PPM from Handheld with Probes

I was able to test the output, which is shown here.

Everything looks good at this point. The next thing to do is connect the voltage and ground pins of the module directly to the batter, they RF out line to the antenna, and the RF good line to ground through a resistor. This is shown below. Note the module is underneath the little breadboard.

TxDuino Prototype

And when running the test program, this is what the airplane does.

And here you can see the waveform generated by the completed TxDuino prototype

TxDuino Prototype Moneyshot

And that's that.

Installing Ubuntu 9.04


These are my install notes for installing Ubuntu 9.04 on a Dell XPS with two NVIDIA 7900 GS cards in SLI and three monitors. Some of this is quite specific to this setup.

I did a no-cd install using the method described here.

To install ubuntu without an install CD, I used the net installer. This was done on a system currently using windows XP, that already had an additional partition for linux.

The first thing I needed to do is download the linux kernel and the net installer. These are in two files linux.bin and initrd.gz. I got them from here.

For future versions, replace "jaunty" with the nickname of the current distribution. After downloading these two files in windows, I placed them in `C:boot. This directory didn't exist yet, so I created it.

The next thing I needed was a boot loader that will allow me to choose which operating system to boot into. I used the open source GRUB for windows. I got it from the sourceforge page. The installer wasn't necessary. Instead I downloaded the grub4dos-0.4.4.zip file from the download page. Then I extracted "gldr" into "C:" and extracted "menu.lst" into "C:bootgrub".

Next, I opened menu.lst with a text editor (notepad++). At the very bottom I added the following

title Ubuntu Installer (hd0,1)
kernel (hd0,1)/boot/linux vga=normal ramdisk_size=14972
root=/dev/rd/1 rw
initrd (hd0,1)/boot/initrd.gz

Note that the numbers above refer to the physical drive and logical partition, so this may change. Check the drive that you're working with.

The next thing I did was I edited boot.ini so that the windows boot loader would give me the option of loading GRUB instead of windows xp. To edit boot.ini run msconfig and click on the edit button for boot.ini. Then I added this line at the end

C:gldr="Start GRUB"

Then I restarted the computer. It stopped after the POST and gave the option of booting into windows or GRUB. I selected GRUB and hit enter. When GRUB loaded I went down to the last optionm "Ubuntu Installer" and hit enter. It then booted into the Ubuntu Net Installer

The installer attempted to use DHCP to connect to the internet. It then asked for a hostname so I put in the static IP address for this machine.

In the last step of the installer it said something about "only the base system is installed" and it asked what other parts to install. I selected something and hit "enter" but that didn't select the item for install, it just advanced to the next step. Thus I ended up with a base installation and just a console prompt. So I installed the rest of the desktop OS with

sudo apt-get install ubuntu-desktop

Then I rebooted to get the desktop. When I was logged in, there was a little green circuit icon on the top right of the desktop. When I clicked this, it prompted me to install the closed-source nvidia drivers. I accepted that and rebooted. I had some problems with the multiple monitors though so I had to edit xorg.conf. Next time, don't reboot, but first open xorg.conf for editing.

sudo gedit /etc/X11/xorg.conf

Then I opened a separate terminal and run

sudo lspci | grep VGA

The first column of what prints out is the "Bus Id" and I used that in xorg.conf need that in a minute. Back in xorg.conf and I edited it as follows.

# nvidia-settings: X configuration file generated by nvidia-settings
# nvidia-settings:  version 1.0  (buildd@palmer)  Sun Feb  1 20:21:04 UTC 2009

Section "ServerLayout"
    Identifier     "Layout0″
    Screen      0  "Screen0 1680 0
    Screen      1  "Screen1″ RightOf "Screen0
    Screen      2  "Screen2″ LeftOf "Screen0
    InputDevice    "Keyboard0″ "CoreKeyboard"
    InputDevice    "Mouse0 "CorePointer"
EndSection

Section "Files"
EndSection

Section "Module"
    Load           "dbe"
    Load           "extmod"
    Load           "type1″
    Load           "freetype"
    Load           "glx"
EndSection

Section "ServerFlags"
    Option         "Xinerama" "1
EndSection

Section "InputDevice"
    # generated from default
    Identifier     "Mouse0″
    Driver         "mouse"
    Option         "Protocol" "auto"
    Option         "Device" "/dev/psaux"
    Option         "Emulate3Buttons" "no"
    Option         "ZAxisMapping" "4 5
EndSection

Section "InputDevice"
    # generated from default
    Identifier     "Keyboard0″
    Driver         "kbd"
EndSection

Section "Monitor"
    # HorizSync source: edid, VertRefresh source: edid
    Identifier     "Monitor2
    VendorName     "Unknown"
    ModelName      "DELL 2007WFP"
    HorizSync       30.0 - 83.0
    VertRefresh     56.0 - 76.0
    Option         "DPMS"
EndSection

Section "Monitor"
    # HorizSync source: edid, VertRefresh source: edid
    Identifier     "Monitor0″
    VendorName     "Unknown"
    ModelName      "DELL 2007WFP"
    HorizSync       30.0 - 83.0
    VertRefresh     56.0 - 76.0
    Option         "DPMS"
EndSection

Section "Monitor"
    # HorizSync source: edid, VertRefresh source: edid
    Identifier     "Monitor1
    VendorName     "Unknown"
    ModelName      "DELL 2007WFP"
    HorizSync       30.0 - 83.0
    VertRefresh     56.0 - 76.0
    Option         "DPMS"
EndSection

Section "Device"
    Identifier     "Device2″
    Driver         "nvidia"
    VendorName     "NVIDIA Corporation"
    BoardName      "GeForce 7900 GS"
    BusID          "PCI:1:0:0
    Screen          1
EndSection

Section "Device"
    Identifier     "Device0″
    Driver         "nvidia"
    VendorName     "NVIDIA Corporation"
    BoardName      "GeForce 7900 GS"
    BusID          "PCI:1:0:0
    Screen          0
EndSection

Section "Device"
    Identifier     "Device1″
    Driver         "nvidia"
    VendorName     "NVIDIA Corporation"
    BoardName      "GeForce 7900 GS"
    BusID          "PCI:6:0:0
EndSection

Section "Screen"
    Identifier     "Screen2″
    Device         "Device2
    Monitor        "Monitor2″
    DefaultDepth    24
    Option         "TwinView" "0
    Option         "metamodes" "DFP-1: nvidia-auto-select +0+0″
    SubSection     "Display"
        Depth       24
    EndSubSection
EndSection

Section "Screen"
    Identifier     "Screen0
    Device         "Device0″
    Monitor        "Monitor0
    DefaultDepth    24
    Option         "TwinView" "0″
    Option         "TwinViewXineramaInfoOrder" "DFP-0
    Option         "metamodes" "DFP-0: nvidia-auto-select +0+0; DFP-0: 1600x1024 +0+0; DFP-0: 1440x900 +0+0; DFP-0: 1400x1050 +0+0; DFP-0: 1360x768 +0+0; DFP-0: 1360x768_60_0 +0+0; DFP-0: 1280x1024 +0+0; DFP-0: 1280x1024_60 +0+0; DFP-0: 1280x960 +0+0; DFP-0: 1152x864 +0+0; DFP-0: 1152x864_75_0 +0+0; DFP-0: 1152x864_70 +0+0; DFP-0: 1152x864_60 +0+0; DFP-0: 1024x768 +0+0; DFP-0: 1024x768_70 +0+0; DFP-0: 1024x768_60 +0+0; DFP-0: 960x600 +0+0; DFP-0: 960x540 +0+0; DFP-0: 896x672 +0+0; DFP-0: 840x525 +0+0; DFP-0: 840x525d70 +0+0; DFP-0: 840x525d60 +0+0; DFP-0: 840x525d60_0 +0+0; DFP-0: 832x624 +0+0; DFP-0: 800x600 +0+0; DFP-0: 800x600d60 +0+0; DFP-0: 800x600_75 +0+0; DFP-0: 800x600_72 +0+0; DFP-0: 800x600_60 +0+0; DFP-0: 800x600_56 +0+0; DFP-0: 800x512 +0+0; DFP-0: 720x450 +0+0; DFP-0: 680x384 +0+0; DFP-0: 680x384d60_0 +0+0; DFP-0: 640x512 +0+0; DFP-0: 640x512d60 +0+0; DFP-0: 640x480 +0+0″
    SubSection     "Display"
        Depth       24
    EndSubSection
EndSection

Section "Screen"
    Identifier     "Screen1
    Device         "Device1″
    Monitor        "Monitor1
    DefaultDepth    24
    Option         "TwinView" "0″
    Option         "metamodes" "nvidia-auto-select +0+0
    SubSection     "Display"
        Depth       24
    EndSubSection
EndSection

Then I needed to install flash in Firefox. I opened a terminal and typed.

sudo apt-get install flash-nonfree

But then I needed to fix a few things because flash is buggy and closed source. I opened a terminal and typed the following.

# Flash also looks for /usr/lib/libesd.so.1
sudo ln -s /usr/lib/libesd.so.0 /usr/lib/libesd.so.1

# Flash expects /tmp/.esd/socket to exist.
sudo mkdir -p /tmp/.esd/
sudo touch /tmp/.esd/socket

And this enabled sound in flash pages.

Equation for the Line of Night on a Map of Earth


Recently a friend of mine expressed an interest in knowin g what the equation was for the night/day line on a rectangular map of the earth. For example:

Line of Night

The line of night is generated by a projection of the Sun's light onto the Earth, which we will consider to be a perfect spheroid with semi-major axis of $latex R_1$ and semi-minor axis of $latex R_2$. If viewed from the "side" the line of night will appear to be a line with slope equal to that of the axial tilt (a.k.a obliquity) of earth. This is the inclination angle of Earth's planetary rotation axis in relation to its orbital plane, usually denoted by $latex \epsilon$

Earth / Sun System

A typical map of the earth is a projection from spherical coordinates onto a cylindrical plane. This is done by drawing a rectangular image where the latitude ($latex \lambda$) is mapped to the x-axis and the longitude ($latex \phi$) is mapped to the y-axis (contrary to the coordinate system used by the digital imaging industry, here we describe the origin as being in the center of the image, with the positive x-axis being to the right, and positive y-axis being up).

Consider then a plane that is perpendicular to the earth's plane of orbit, and perpendicular to the vector pointing from the center of the earth toward the sun, and that intersects the earth at it's center. The line of night is the intersection of this plane with the Earth. From this observation we can develop an equation for the line of night in terms of $latex \phi$ and $latex \lambda$.

Lattitude and Longitude

[Line of Night

If we place a coordinate system on the surface of the earth, centered at the intersection of the line of night with the equator with the positive x-axis corresponding to positive $latex \lambda$ and positive y-axis corresponding to positive $latex \phi$ we can describe the line-of-night plane as

$$y = s x$$

where the slope $latex s$ is related to the obliquity as

$$ s = \frac{1}{\tan(\epsilon)} $$

where $latex \epsilon$ is the inclination of the axial tilt of the Earth's axis at solstice ($latex \epsilon_0$) scaled according to the position of the Earth in it's orbit

$$ \epsilon = \epsilon_0 \cos( \frac{day}{360} ) $$

The inclination of the earth's axis at solstice is $latex \epsilon_0 = 23.5^\circ $

From the parameterization of a spheroid, we know that

$$ x = R_1 \sin( \lambda ) \cos( \phi ) $$

$$ y = R_2 \sin( \phi ) $$

Combining these with the equation above for the line-of-night plane we see that

$$ y = \frac{1}{\tan( \epsilon )} x $$

$$ R_2 \sin( \phi ) = \frac{1}{\tan( \epsilon )} R_1 \sin( \lambda ) \cos( \phi ) $$

$$ \tan( \phi ) = \frac{1}{\tan( \epsilon )} \frac{R_1}{R_2} \sin( \lambda ) $$

$$ \phi = \text{tan}^-1 \left( \frac{1}{\tan( \epsilon )} \frac{R_1}{R_2} \sin( \lambda ) \right) $$

Plotting the value of the longitude against latitude over a range of axial inclinations from $latex \epsilon = 0$ to $latex \epsilon = 23.5^\circ$, and assuming $latex R_1 = R_2$ results in the following:

Plot of Resulting Line-of-Night Equation

Which has exactly the shape we were looking for.

Point on an Ellipsoid of Minimum Distance to Another Point in Space


This is a response to a question posted on the irrlicht forums at here

PDF file:
Closest point on an ellipsoid to a point in space (to come)

If one wants to know the closest point on an ellipsoid to another point in space, it can be found using constrained nonlinear optimization. It was stated in the forum thread that there is no closed-form solution to this problem, though I find that rather surprising. I will update this post if I can confirm (or confirm an opposition) to that fact.

In any case, given a point $latex \vec{p}$ in space, where

$$ \vec{p} = \begin{bmatrix} p_x && p_y && p_z \end{bmatrix} $$

we want to find a point $latex \vec{q}$ that minimizes the euclidean distance to $latex \vec{p}$ subject to the contraint that it is on an ellipsoid with the
following functional description

$$ \frac{x^2}{a^2} + \frac{y^2}{b^2} + \frac{z^2}{c^2} = 1 $$

This corresponds to the minimization of the following cost function (the "cost" is the distance to the target point)

$$ J(\vec{q}) = \sqrt{ (p_x-q_x)^2 + (p_y-q_y)^2 + (p_z-q_z)^2 } $$

subject to the following constraint (which just says that it must be on the ellipsoid)

$$ f(\vec{q}) = \frac{q_x^2}{a} + \frac{q_y^2}{b} + \frac{q_z^2}{c} - 1 = 0 $$

Now there are many ways to solve this problem. I will outline two of them here known as steepest descent (or gradient search) and the leGrange method.

Steepest Descent

First comes steepest decent, because it is the easiest to understand graphically. The goal here is that, starting at some point on the ellipse, we will look at how the value of the cost function changes depending on the direction we move, and choose to move in the direction that most improves the cost. To begin with, let us redefine the problem in terms of two independent variables, and let the third be dependent. In the following formulation, I will choose $latex q_z$ to be the dependent variable, but you should note that this is not always the best choice. A good choice of dependent variable is to look at the coordinates of the target point and pick the one that is furthest from the origin (to avoid numerically ill-conditioned situations).

Using $latex q_z$ as the dependent variable we can reformulate the constraint as the following

$$q_z = \pm c \sqrt{ 1 - \frac{q_x^2}{a^2} - \frac{q_y^2}{b^2} }$$

Given this, we exam the effect that changing each of our independent variables has on our cost function, by calculating the partial derivatives. First though, we note that the minimum of a square root expression will occur at the same point as the minimum of the argument, so we will throw away the square root part and use

$$ J = (p_x-q_x)^2 + (p_y-q_y)^2 + (p_z-q_z)^2 $$

for the derivatives instead.

$$ \frac{ \partial J }{ \partial q_x } = -2(p_x - q_x) + -2(p_z - q_z)\frac{ \partial q_z }{ \partial q_x }$$
$$ \frac{ \partial J }{ \partial q_y } = -2(p_y - q_y) + -2(p_z - q_z)\frac{ \partial q_z }{ \partial q_y }$$

And given the equation above for $latex q_z$ we can determine the following

$$ \frac{ \partial q_z }{ \partial q_x } = \pm c \frac{1}{2} \sqrt{ 1 - \frac{q_x^2}{a^2} - \frac{q_y^2}{b^2} }^{-\frac{1}{2}} \left( -2 \frac{q_x}{a^2}\right)$$

which simplifies to

$$ \frac{ \partial q_z }{ \partial q_x } = \mp \frac{c}{a^2} \frac{q_x}{q_z}$$

and similarly

$$ \frac{ \partial q_z }{ \partial q_y } = \mp \frac{c}{b^2} \frac{q_y}{q_z}$$

Note that the sign of these derivatives is opposite the sign of $latex q_z$. Therefore we see the following

$$ \frac{ \partial J }{ \partial q_x } = -2(p_x - q_x) + 2 sign(q_z) (p_z - q_z) \left( \frac{c}{a^2} \frac{q_x}{q_z} \right)$$
$$ \frac{ \partial J }{ \partial q_y } = -2(p_y - q_y) + 2 sign(q_z) (p_z - q_z) \left( \frac{c}{b^2} \frac{q_y}{q_z} \right)$$

where $latex sign(x)$ is the signum function and returns $latex -1$ if the argument is less than zero ($latex x < 0$) or $latex 1$ if the argument is greater than zero ($latex x > 0$).

Now we solve for the minimum point using the following algorithm

  1. start at some $latex q_x$ and $latex q_y$
  2. calculate $latex \pm q_z$
  3. pick whichever of the two is closer by picking the one with the same sign as $latex p_z$
  4. calculate the cost $latex J$
  5. if $latex J$ of the current iteration is within some tolerance of $latex J$ at the previous iteration, stop here
  6. calcluate $latex \frac{\partial J}{\partial q_x}$ and $latex \frac{\partial J}{\partial q_y}$
  7. calculate $latex \Delta x$ and $latex \Delta y$ from $latex J$, $latex \frac{\partial J}{\partial q_x}$, and $latex \frac{\partial J}{\partial q_y}$
  8. update $latex q_x$ and $latex q_y$ using $latex \Delta x$ and $latex \Delta y$
  9. return to (1)

Here is a matlab code that demonstrates this algorithm. Note that in this code, for the initial guess I use the intersection of the line segment between the origin of the ellipse and $latex \vec{p}$ with the surface of the ellipse.

    clc;
    clear;

    % define the ellipsoid parameters
    a   = 3;
    b   = 7;
    c   = 9;

    % createa a 100 x 100 point mesh for displaying the ellipsoid
    x   = linspace( -a, a, 50 );
    y   = linspace( -b, b, 50 );

    % calculate the z coordinates
    [X,Y]   = meshgrid(x,y);
    Zp      = c*sqrt(1 - X.*X/(a*a) - Y.*Y/(b*b) );
    Zm      = -Zp;

    % create the point we're going to calculate the closest too
    px = 3;
    py = 4;
    pz = 8;


    % create an initial guess by finding the intersection of the ray from the
    % centroid to point with the ellipsoid surface
    k = sqrt( 1 / ( (px^2)/(a^2) + (py^2)/(b^2) + (pz^2)/(c^2) ) );

    qx = k*px;
    qy = k*py;
    qz = k*pz;

    % this is the scale factor for our movement along the ellipse, this can
    % start off quite large because it will be reduced as needed
    s = 1;


    % this is the main loop, you'll want to put some kind of tolerance based
    % terminal condition, i.e. when the cost function (distance) is only 
    % bouncing back and forth between a tolerance or something more tailored to
    % your application, I will stop when the distance is +/- 1%
    i = 1;
    while 1
        % calculate the z for our given x and y
        qz_plus     = c * sqrt( 1 - qx*qx/(a*a) - qy*qy/(b*b) );
        qz_minus    = - qz_plus;

        % we want the one that get's us closest to the target so 
        if( abs(pz - qz_plus) < abs(pz - qz_minus) )
            qz = qz_plus;
        else
            qz = qz_minus;
        end



        % calculate the current value of the cost function
        J           = sqrt( (px - qx)^2 + (py - qy)^2 + (pz - qz)^2 );

        % store the current values for the plots
        qxplot(i)   = qx;
        qyplot(i)   = qy;
        qzplot(i)   = qz;
        Jplot(i)    = J;

        % check to see if we overshot the goal or jumped off the ellipsoid
        if( i > 1 )
            % if we jumped off the ellipsoid or overshot the minimal cost
            if( imag(qz) ~= 0 || J > Jplot(i-1) );
                % then go back to the previous position and use a finer 
                % step size
                qx = qxplot(i-1);
                qy = qyplot(i-1);
                qz = qzplot(i-1);
                J  = Jplot(i-1);
                s = s/10;
                i = i-1;

                % if we did just jump over the actual minimum, let's check to 
                % see how confident we are at this point, if we're with in 
                % 1% there's no need to continue
                if( i > 3 )
                    if( abs( (Jplot(i-1) - Jplot(i-3)) / Jplot(i-3) ) < 0.001 )
                        break;
                    end
                end
            end
        end

        % calculate the gradient of the cost with respect to our control
        % variables; since we divide by qz in the second term we first need
        % to check whether or not qz is too close to zero; if it is, we know
        % that the second term should be zero
        if( qz < 1e-10 )
            dJdx = -2*(px - qx);
            dJdy = -2*(py - qy);
        else
            dJdx    = -2*(px - qx) + 2*sign(qz)*(pz - qz)*( c/(a*a) * qx/qz );
            dJdy    = -2*(py - qy) + 2*sign(qz)*(pz - qz)*( c/(b*b) * qy/qz );
        end

        % calculate the update \vector that we will move along
        dx          = -J/dJdx;
        dy          = -J/dJdy;

        % calculate the magnitude of that update \vector
        magnitude   = sqrt( dx*dx + dy*dy );

        % normalize our update \vector so that we don't shoot off at an
        % uncontrolable rate
        dx          = s * (1/magnitude) *dx;
        dy          = s * (1/magnitude) *dy;

        % update the current position 
        qx          = qx + dx;
        qy          = qy + dy;

        % increment the index 
        i = i+1;
    end;

    index = 1:i+1;

    % generate a little report 
    message = sprintf( ...
        ['The closest point was found at [%f, %f, %f]', ...
         'with a distance of %f and a confidence of +/- %f %%'], ...
        qx, qy, qz, J, 100*abs((Jplot(i+1)-Jplot(i-1)) / Jplot(i-1)) );
    disp(message);

    % now we'll verify that the one we found is the closest
    Jmap_plus   = sqrt( (px - X).*(px - X) + (py - Y).*(py - Y) + (pz - Zp).*(pz - Zp ) );
    Jmap_minus  = sqrt( (px - X).*(px - X) + (py - Y).*(py - Y) + (pz - Zp).*(pz - Zm ) );

    [Jmin_plus_row, imin_plus_row] = min( Jmap_plus );
    [Jmin_plus_col, imin_plus_col] = min( Jmin_plus_row );

    imin = imin_plus_row(imin_plus_col);
    jmin = imin_plus_col;

    xmin    = X(imin, jmin);
    ymin    = Y(imin, jmin);
    zmin    = Zp(imin, jmin);
    Jmin    = Jmin_plus_col;
    Jmin2   = Jmap_plus(imin,jmin);

    message = sprintf( ... 
        ['The closest point by searching the mesh was at [%f, %f, %f]' ...
         'and has a distance of %f (%f)'], xmin, ymin, zmin, Jmin, Jmin2 );
    disp(message);


    % and draw some pretty pictures

    close(figure(1));
    figure(1);
    grid on;
    hold on;
    mesh(X,Y,real(Zp), 'FaceAlpha', 0.5);
    mesh(X,Y,real(Zm), 'FaceAlpha', 0.5);
    plot3(px,py,pz,'ko');
    plot3(qxplot,qyplot,qzplot,'ko');
    plot3(qxplot,qyplot,qzplot,'k-');
    plot3([px qx], [py qy], [pz qz], 'b-');
    xlabel('X');
    ylabel('Y');
    zlabel('Z');
    title('steepest decent method');
    hold off;

    close(figure(2));
    figure(2);
    plot(index, Jplot, 'g-');
    xlabel( 'iteration' );
    ylabel( 'distance (cost)' );

    close(figure(3));
    figure(3);
    mesh(X,Y,real(Jmap_plus), 'FaceAlpha', 0.5);
    xlabel('X');
    ylabel('Y');
    zlabel('Z');
    title('distance to target');

The output of this script is

The closest point was found at [1.321810, 3.187960, 6.962411]with a distance
of 2.133617 and a confidence of +/- 0.000000 %  
The closest point by searching the mesh was at [1.285714, 3.285714,
6.948103]and has a distance of 2.134354 (2.134354)  

You can see the work done by the gradient search in the first figure, which is shown below:

Result of Steepest Descent Algorithm

Full Size

Detail of Steepest Descent

Full Size

I know this doesn't look like the closest point, but, as you can see, the script actually searches the mesh for the lowest value point and comes up with essentially the same point. To demonstrate, here is a 3d plot of the distance of each point on the upper half of the ellipse to the target point, with the minimum point labeled. This is the output of the script to figure 3. The reason it doesn't look right is because of how matlab squashes everything to fit the window.

[Distance of Point on Ellipse to Target Point

Full Size

And finally here is a plot of the distance calculated at each iteration of the loop. As you can see, it doesn't take many to find a value that is pretty damn good.

ellipse_steepest_distance

Full Size

Steepest decent works pretty well right? Well, not always. The initial guess used in the code above will save your ass, but steepest descent can suffer from problems of local minima. Look what happens when I try an initial guess on the other side of the ellipsoid.

Failure of Steepest Descent Algorithm

Full Size

The algorithm get's caught in a local minimum, and converges to a very precise, very incorrect solution.

LeGrange Method

This method is a little more elegant, but still ends up with a search algorithm. However, in this case we're only doing a one dimensional search for root-finding of a polynomial, which is a quite well documented problem.

In this method, we note that, since the constraint must be satisfied everywhere, if we know the minimal cost $latex J^$ is achieved at a point $latex \vec{q}^$, then we know that

$$ J^ = \sqrt{ (p_x-q^_x)^2 + (p_y-q^_y)^2 + (p_z-q^_z)^2 } $$

and we also know that if the augmented cost function is defined as

$$ L = J(\vec{q}) + \lambda f(\vec{q}) $$

which expands to

$$ L = \sqrt{ (p_x-q_x)^2 + (p_y-q_y)^2 + (p_z-q_z)^2 } + \lambda \left[ \frac{q_x^2}{a^2} + \frac{q_y^2}{b^2} + \frac{q_z^2}{c^2} - 1 \right] $$

then $latex L^ = J^$ and is also at $latex \vec{q}^*$ since the latter term must be zero if the constraint is to be satisfied. We can then minimize the cost function as we usually do, by taking the derivative with respect to the independent variables, setting them to zero, and solving the system for the unknowns.

First though, we note that the minimum of a square root expression will occur at the same point as the minimum of the argument, so we will throw away the square root part and solve instead for the minimum of

$$ L = (p_x-q_x)^2 + (p_y-q_y)^2 + (p_z-q_z)^2 + \lambda \left[ \frac{q_x^2}{a} + \frac{q_y^2}{b} + \frac{q_z^2}{c} - 1 \right] $$

We take the partial derivative of the augmented cost with respect to
the four unknowns and see that

$$ \frac{\partial{L}}{\partial{q_x}} = -2(p_x - q_x) + \frac{2 \lambda q_x}{a^2} $$
.
$$ \frac{\partial{L}}{\partial{q_y}} = -2(p_y - q_y) + \frac{2 \lambda q_y}{b^2} $$
.
$$ \frac{\partial{L}}{\partial{q_z}} = -2(p_z - q_z) + \frac{2 \lambda q_z}{c^2} $$
.
$$ \frac{\partial{L}}{\partial{\lambda}} = \frac{q_x^2}{a^2} + \frac{q_y^2}{b^2} + \frac{q_z^2}{c^2} - 1 $$

Setting these all to zero we can solve for the following

$$ q_x = \frac{p_x}{\lambda + \frac{1}{a^2}} = \frac{a^2 p_x}{\lambda + 1} = 0 $$
.
$$ q_y = \frac{p_y}{\lambda + \frac{1}{b^2}} = \frac{b^2 p_y}{\lambda + 1} = 0 $$
.
$$ q_z = \frac{p_z}{\lambda + \frac{1}{c^2}} = \frac{c^2 p_z}{\lambda + 1} = 0 $$

Substituting back into the fourth of the differentials we get that

$$ \frac{\partial{L}}{\partial{\lambda}} = \frac{a^2 p_x^2}{ \left(\lambda + a^2 \right)^2} + \frac{b^2 p_y^2}{ \left(\lambda + b^2 \right)^2} + \frac{c^2 p_z^2}{ \left(\lambda + c^2 \right)^2} - 1 = 0 $$

which simplifies to

$$ \frac{\partial{L}}{\partial{\lambda}} = \begin{matrix} & a^2 p_x^2\left( \lambda + b^2 \right)^2 \left( \lambda + c^2 \right)^2 \ & + \ b^2 p_y^2 \left( \lambda + a^2 \right)^2 \left( \lambda + c^2 \right)^2 \ & + \ c^2 p_z^2 \left( \lambda + a^2 \right)^2 \left( \lambda + b^2 \right)^2 \ & - \left( \lambda + a^2 \right)^2 \left( \lambda + b^2 \right)^2 \left( \lambda + c^2 \right)^2 = 0 \end{matrix} $$

Note that this is a sixth order polynomial in $latex \lambda$. Use your favorite zero-finding algorithm (same as a root finding algorithm) to solve for the zeros of this polynomial, (of which there should only be two as there can only be one closest point, and one furthest point on an ellipsoid unless the target point is along one of the axes), and you will find the minimum and the maximum of your cost function. Once the solutions of $latex \lambda$ are known, plug them into the equations above to solve for $latex q_x$, $latex q_y$, and $latex q_z$. Then plug those into the cost function, and pick the smaller of the two.

For the example problem used in the steepest descent code above, the polynomial looks like this:

[Polynomial For Legrange Multiplier

Full Size

Here is a detail view:

Polynomial for LeGrange Multiplier Detail

Full Size

As you can see, there are two solutions. Below is a matlab code that generates the above graphs, and solves for the minimum point using a binary search. Note that the search assumes that one of the values of \lambda that makes the polynomial zero is to the left of the origin, and the other is to the right. That is true for this example buy may not be in general (I didn't take the time to work out why that may be the case ) so you may have to modify the search to look for both roots on both sides of the origin.

    clc;
    clear;

    % define the ellipsoid parameters
    a   = 3;
    b   = 7;
    c   = 9;

    % createa a 100 x 100 point mesh for displaying the ellipsoid
    x   = linspace( -a, a, 50 );
    y   = linspace( -b, b, 50 );

    % calculate the z coordinates
    [X,Y]   = meshgrid(x,y);
    Zp      = c*sqrt(1 - X.*X/(a*a) - Y.*Y/(b*b) );
    Zm      = -Zp;

    % create the point we're going to calculate the closest too
    px = 3;
    py = 4;
    pz = 8;

    \lambda = -160:0.01:20;
    a2 = a*a;
    b2 = b*b;
    c2 = c*c;

    px2 = px*px;
    py2 = py*py;
    pz2 = pz*pz;

    lpa = \lambda + a*a;
    lpb = \lambda + b*b;
    lpc = \lambda + c*c;

    P = a2*px2*lpb.*lpb.*lpc.*lpc +b2*py2*lpa.*lpa.*lpc.*lpc +c2*pz2*lpa.*lpa.*lpb.*lpb -lpa.*lpa.*lpb.*lpb.*lpc.*lpc;

    x = a2*px./lpa;
    y = b2*py./lpb;
    z = c2*pz./lpc;

    J = sqrt( (px-x).*(px-x) +  (py-y).*(py-y) +  (pz-z).*(pz-z) );

    close(figure(1));
    figure(1);
    plot( \lambda, P, 'b-', 'linewidth', 2 );
    xlabel('$\lambda$', 'Interpreter', 'LaTex');
    ylabel('$\frac{\partial L}{ \partial \lambda }$', 'Interpreter', 'LaTex');
    title('value of the polynomial');

    close(figure(2));
    figure(2);
    plot( \lambda, P, 'b-', 'linewidth', 2 );
    axis([-160,20,0, 1e10]);
    xlabel('$\lambda$', 'Interpreter', 'LaTex');
    ylabel('$\frac{\partial L}{ \partial \lambda }$', 'Interpreter', 'LaTex');
    title('detail of value of the polynomial');


    % start search by looking at the value of the polynomial at zero
    \lambda = 0;

    lpa = \lambda + a*a;
    lpb = \lambda + b*b;
    lpc = \lambda + c*c;

    P = a2*px2*lpb.*lpb.*lpc.*lpc +b2*py2*lpa.*lpa.*lpc.*lpc +c2*pz2*lpa.*lpa.*lpb.*lpb -lpa.*lpa.*lpb.*lpb.*lpc.*lpc;

    % we need to store the start sign so that we know what sign we're looking
    % for in order to determine a zero crossing
    startSign = sign(P);

    % start our search in one direction
    i = 1;

    d\lambda = 1;
    lhs\lambda = 0;
    % start the forward search, we're looking for the first value we can 
    % find that has a sign opposite that at the zero point
    while d\lambda < realmax/10
        rhs\lambda   = d\lambda;
        lpa         = rhs\lambda + a*a;
        lpb         = rhs\lambda + b*b;
        lpc         = rhs\lambda + c*c;

        P       =   a2*px2*lpb.*lpb.*lpc.*lpc + ...
                    b2*py2*lpa.*lpa.*lpc.*lpc + ...
                    c2*pz2*lpa.*lpa.*lpb.*lpb - ...
                    lpa.*lpa.*lpb.*lpb.*lpc.*lpc;

        % if we've discovered a sign change we can stop searching
        if( sign(P) ~= startSign )
            lhs\lambda = d\lambda/10;
            break;

        % if not, we need to grow the increment
        else
            d\lambda = d\lambda * 10;
        end
    end


    % now we can start a bisection search using the lhs and rhs that we've
    % just determined, which are exactly one order of magnitude apart
    error = (rhs\lambda - lhs\lambda)/lhs\lambda;
    while( abs(error) > 1e-2)
        \lambda = (lhs\lambda + rhs\lambda)/2;
        lpa         = \lambda + a*a;
        lpb         = \lambda + b*b;
        lpc         = \lambda + c*c;

        P       =   a2*px2*lpb.*lpb.*lpc.*lpc +b2*py2*lpa.*lpa.*lpc.*lpc +c2*pz2*lpa.*lpa.*lpb.*lpb -lpa.*lpa.*lpb.*lpb.*lpc.*lpc;

        if( sign(P) ~= startSign )
            rhs\lambda = \lambda;
        else
            lhs\lambda = \lambda;
        end

        error = (lhs\lambda - rhs\lambda)/lhs\lambda;
    end

    % store the found value
    \lambda1 = (lhs\lambda + rhs\lambda)/2;

    % now we search in the other direction;
    d\lambda = 1;
    rhs\lambda = 0;
    % start the forward search, we're looking for the first value we can 
    % find that has a sign opposite that at the zero point
    while d\lambda < realmax/10
        lhs\lambda   = -d\lambda;
        lpa         = lhs\lambda + a*a;
        lpb         = lhs\lambda + b*b;
        lpc         = lhs\lambda + c*c;

        P       =   a2*px2*lpb.*lpb.*lpc.*lpc + ...
                    b2*py2*lpa.*lpa.*lpc.*lpc + ...
                    c2*pz2*lpa.*lpa.*lpb.*lpb - ...
                    lpa.*lpa.*lpb.*lpb.*lpc.*lpc;

        % if we've discovered a sign change we can stop searching
        if( sign(P) ~= startSign )
            rhs\lambda = -d\lambda/10;
            break;

        % if not, we need to grow the increment
        else
            d\lambda = d\lambda * 10;
        end
    end

    % now we can start a bisection search using the lhs and rhs that we've 
    % just determined, which are exactly one order of magnitude apart
    error = (rhs\lambda - lhs\lambda)/lhs\lambda;
    while( abs(error) > 1e-2)
        \lambda = (lhs\lambda + rhs\lambda)/2;
        lpa         = \lambda + a*a;
        lpb         = \lambda + b*b;
        lpc         = \lambda + c*c;

        P       =   a2*px2*lpb.*lpb.*lpc.*lpc +b2*py2*lpa.*lpa.*lpc.*lpc +c2*pz2*lpa.*lpa.*lpb.*lpb -lpa.*lpa.*lpb.*lpb.*lpc.*lpc;

        if( sign(P) ~= startSign )
            lhs\lambda = \lambda;
        else
            rhs\lambda = \lambda;
        end

        error = (lhs\lambda - rhs\lambda)/lhs\lambda;
    end

    % store the found value
    \lambda2 = (lhs\lambda + rhs\lambda)/2;

    x1 = a2*px/(\lambda1+a2);
    y1 = b2*py/(\lambda1+b2);
    z1 = c2*pz/(\lambda1+c2);

    x2 = a2*px/(\lambda2+a2);
    y2 = b2*py/(\lambda2+b2);
    z2 = c2*pz/(\lambda2+c2);

    J1 = sqrt( (px-x1)*(px-x1) + (py-y1)*(py-y1) + (py-y1)*(py-y1) );
    J2 = sqrt( (px-x2)*(px-x2) + (py-y2)*(py-y2) + (py-y2)*(py-y2) );

    % print a little report
    message = sprintf( ['found zero crossings at \lambda = %f and %f', ...
                        'which corresponds to the points [%f, %f, %f]','and [%f, %f, %f] with costs of %f and %f'],\lambda1, \lambda2, x1, y1, z1, x2, y2, z2, J1, J2 );
    disp(message);

    if( J1 < J2 )
        message = sprintf( 'the point of minimum distance is [%f, %f, %f]', ...
                            x1, y1, z1 );
    else
        message = sprintf( 'the point of minimum distance is [%f, %f, %f]', ...
                            x2, y2, z2 );
    end

    disp(message);

The output of this script is

found zero crossings at \lambda = 11.801758 and -155.810547which corresponds
to the points [1.297967, 3.223591, 6.982626]and [-0.183910, -1.835025,
-8.661880] with costs of 2.025472 and 8.844903  
the point of minimum distance is [1.297967, 3.223591, 6.982626]  

Which matches with the previous method. Sweet.

As a brief note, if the target point is inside the ellipse and along one of it's axes, there may be an infinite number of solutions (all points on a ring around the ellipse located at the coordinate of the target point will be the same distance away). If this is a possibility in your application, make sure you watch out for it and decide what the "right" solution is.

I hope this helps. Sorry I didn't write the code in c++ but it took long enough to generate the results without doing that. Hopefully you can figure out from the matlab code exactly what you need to. Matlab script isn't too different from C++.

-Cheshirekow

Instanciate Objects of Unknown Type from Their Parent Interface


This is based on my previous two posts on Static Interfaces in C++ and Keep Track of and Enumerate All Sub-classes of a Particular Interface. The idea is that I want my code to be extensible in the feature without requiring any re-writing of the current code base. The code base operates on generic objects via their interfaces, so as long as newly-coded classes properly extend those interfaces, the program should know how to handle them. The problem is, how can we write the program in such a manner that a user interface can enumerate available options for implementations of a particular interface, and how can we instantiate those objects?

In Keep Track of and Enumerate All Sub-classes of a Particular Interface I showed how to maintain a registry of classes deriving from a given interface, which handles the first problem, but there is a limitation in that all of these classes must provide a factory method that takes no parameters (void input). I decided that, for my project, this was not acceptable and I needed a way to define the creation parameters as part of the factory methods, whereas the creation parameters may be different for particular interfaces.

In Static Interfaces in C++ I showed how we can enforce the requirement of a static method in derived classes with a particular signature using a template interface.

In this post I will combine the two so that we can create a registry of classes that inherit from a particular interface, and provide a static factory method for creating objects of that interface, using a particular creation method signature unique to that interface. The registry will pair class names with function pointers that match the specific signature of the interface the class is being registered for.

Disclaimer: I do not claim this is the "best" way to handle this issue. This is just what I came up with. It happens to be pretty involved and overly indirect, which means it's probably bad design. It is, however, an extremely interesting exercise in generic programming.

Prequil: the code will require these later so there they are:

    /**
     *  file          RegistryTest.cpp
     *  date:      Aug 14, 2009
     *  brief:
     *
     *  detail:
     */


    #include 
    #include 
    #include 
    #include 

    using namespace std;

Ok, so lets begin. First let's define a couple of interfaces that we're interested in.

    class InterfaceA{};
    class InterfaceB{};
    class InterfaceC{};

Now we create a template class whose sole purpose is to create a per-interface typedef of the function signature that is necessary for instantiating and object of that class. Is it really possible that all sub-objects can be instantiated with the same parameters? If that's the case, shouldn't they all be combined into a single class that just contains that information as private members? Probably, but in my case these parameters are more like a "bare minimum" for instantiation, and then many more parameters are set by the user. It makes sense to me, I promise. If it doesn't to you, you don't have to use this.

    template< typename InterfaceType >
    class Factory
    {
        public:
            typedef InterfaceType*(*Creator)(void);
    };

Creator is now a typedef that aliases a function pointer that takes no parameters. Wait, isn't that what we had before? Yes, but now we make a couple of template specializations to define the different signatures for our specific interfaces. These specializations would normally be in the file that contained the interface declaration.

    /// specializations can define other creators, this one requires an int
    template<> class Factory
    {
        public:
            typedef InterfaceB*(*Creator)(int);
    };

    /// specializations can define other creators, this one requires an int, a
    /// bool, and a char
    template<> class Factory
    {
        public:
            typedef InterfaceC*(*Creator)(int,bool,char);
    };

Cool. Now we create a static interface that enforces it's derivative classes to contain a static method called createNew which can be used to instantiate a new object of that interface. We can use the typedef we just created to make the function signature generic for this template (or specific to individual instantiations of it).

    template
    class IStaticFactory
    {
        public:
            IStaticFactory()
            {
                typename Factory::Creator check = ClassType::createNew;
                check = check;
            }
    };

Still following? Good. Now we define the registry class template, which maps the class name of a derived class to a function pointer with an interface- specific signature that serves as a static factory for objects of the derived class, returning a pointer to that object of the type of the interface. See my previous post for details on this class.

    template 
    class Registry
    {
        private:
            std::map< std::string, typename Factory::Creator > m_creatorMap;


            Registry(){}

        public:

            static Registry& getInstance();


            bool registerClass( const std::string& name,
                                 typename Factory::Creator creator );


            std::set  getClassNames();


            typename
            Factory::Creator
            Registry::getCreator(
                    std::string className );
    };



    // A convient macro to compact the registration of a class
    #define RegisterWithInterface( CLASS, INTERFACE )                   
    namespace                                                           
    {                                                                   
        bool dummy_ ## CLASS =                                          
        Registry::getInstance().registerClass(     
                #CLASS, CLASS::createNew );                               
    }


    template 
    Registry& Registry::getInstance()
    {
        static Registry    registry;
        return registry;
    }


    template 
    bool Registry::registerClass( const std::string& name,
                            typename Factory::Creator creator )
    {
        m_creatorMap[name] = creator;
        return true;
    }

    template 
    std::set  Registry::getClassNames()
    {
        std::set   keys;

        typename
        std::map< std::string, InterfaceType* (*)(void) >::iterator pair;

        for( pair = m_creatorMap.begin(); pair != m_creatorMap.end(); pair++)
            keys.insert( pair->first );

        return keys;
    }

    template 
    typename
    Factory::Creator
    Registry::getCreator(
            std::string className )
    {
        return m_creatorMap[className];
    }

The difference between this and the Registry in my previous post, is that this time the registry uses the generic Factory<InterfaceType>::Creator typedef to define the function pointer. This way, that pointer is forced to have the specific signature. Sweet!

Now lets write some derived classes of those interfaces.

    class   DerivedA : public InterfaceA,
                        public IStaticFactory
    {
        public:
            static InterfaceA*  createNew(){ return (InterfaceA*)1; }
    };


    RegisterWithInterface(DerivedA, InterfaceA);



    class   DerivedB : public InterfaceB,
                        public IStaticFactory
    {
        public:
            static InterfaceB*  createNew(int a){ return (InterfaceB*)2; }
    };


    RegisterWithInterface(DerivedB, InterfaceB);



    class   DerivedC : public InterfaceC,
                        public IStaticFactory
    {
        public:
            static InterfaceC*  createNew(int a, bool b, char c){ return (InterfaceC*)3; }
    };

    RegisterWithInterface(DerivedC, InterfaceC);

These classes are basically dummies, but inheriting from IStaticFactory… the compiler will enforce that they contain the static method createNew with the proper signature. Notice that InterfaceA uses the default template so the static factory in DerivedA takes no parameters, while InterfaceB and InterfaceC have specializations so the static factories in DerivedB and DerivedC have their respective parameters. Since this is just an example, the methods don't actually create new objects they just return pointers, but in reality this is where we would use new DerivedA(…) and so on.

Well that's it. Pretty cool huh? The compiler will enforce all this stuff for us so we can actually say to ourselves when we write new implementations months from now "If it compiles, it will be compatible."

Lastly, here's a little test case to run

    int main()
    {
        DerivedA    a;
        DerivedB    b;
        DerivedC    c;

        InterfaceA* pA;
        InterfaceB* pB;
        InterfaceC* pC;

        Factory::Creator makesObjectOfA =
                Registry::getInstance().getCreator("DerivedA");
        pA = (*makesObjectOfA)();

        Factory::Creator makesObjectOfB =
                Registry::getInstance().getCreator("DerivedB");
        pB = (*makesObjectOfB)(1);

        Factory::Creator makesObjectOfC =
                Registry::getInstance().getCreator("DerivedC");
        pC = (*makesObjectOfC)(1,false,'a');

        cout << "pA: " << pA << "n";
        cout << "pB: " << pB << "n";
        cout << "pC: " << pC << "n";

        return 0;
    }

Output of the Test Program

Static Interfaces in C++


I remember looking around a few weeks ago for how to make a "static interface" in c++. Basically, I wanted a way to use the compiler to enforce that a class had certain static functions. Almost all of the internet resources I found basically said "Why would you ever want to do that; you don't really want to do that; you probably have bad design" and so on… continuously begging the question. Of course, they were right: the design was bad and that wasn't really what I wanted to do. Well, never the less, I still managed to think of a way to create a sort of static interface using a template class.

The strategy is to define a template class that uses the static methods of the template parameter class. That way, as long as the template is instantiated, the compiler will complain unless we have provided those static functions. We can ensure that the template is instantiated and enforce the inheritance idea by making the derived class extend from the template class we wrote to enforce those static methods.

Here is an example. We can create the static interface by declaring a class template that uses the functions we want to enforce as part of the interface.

    template < typename T >
    class StaticInterface
    {
        public:
            StaticInterface()
            {
                int(*fooCheck)(int)     = T::foo;
                bool(*barCheck)(bool)   = T::bar;
            }
    };

By assigning T::foo and T::bar to function pointers, we are saying, implicitly, that whatever class is provided as a parameter to this template must have a static method called foo and a static method called bar and, furthermore, that those static methods must have
the same signature as the function pointers we stuff them into.

By putting this code inside the constructor of the class, we know that this method of the template will be instantiated, even if we don't explicitly use it later in the code, as long as we derive from this class somewhere. So then, the last question is, where can we derive from it?

Well, in the class that we want to inherit the interface of course!

    class DerivedClass : public StaticInterface
    {
        public:
            static int foo(int  param){ return 10; }
            static bool bar(bool param){ return 20; }
    };

The DerivedClass constructor implicitly calls the StaticInterface constructor, which assigns the function pointers fooCheck and barCheck to the address of the functions DerivedClass::foo and DerivedClass::bar. As a result, if we forget the bar function in the DerivedClass the compiler will choke with an error. g++ says the following:

src/poc/test/StaticInterfaceTest.cpp: In constructor
`StaticInterface::StaticInterface() [with T = DerivedClass]':  
src/poc/test/StaticInterfaceTest.cpp:41: instantiated from here  
src/poc/test/StaticInterfaceTest.cpp:20: error: `bar' is not a member of
`DerivedClass'

Pretty cool huh?

As a final note, please consider this "an interesting observation" and not necessarily a "great design choice". As I said, I decided against actually trying to utilize this idea in my project, and I urge you think carefully before about yours before trying to use it yourself.

Keep Track of and Enumerate All Sub-classes of a Particular Interface


Here's another interesting C++ problem I encountered today. As is often the case, I found my answer by asking over at Stack Overflow. The question can be found here: http://stackoverflow.com/questions/1260954/how-can-i-keep-track- of-enumerate-all-classes-that-implement-an-interface

The issue was that I was writing a simulation program where I knew I would eventually want to simulate multiple different vehicles, with multiple different controllers, and multiple different estimators. Naturally, this led me to define an interface for each of these things, but the problem was that I only really want to implement a few subclasses of these interfaces now. In particular, I only have two vehicles, 2 controllers, and 1 estimator I'm interested in completing now, but I will probably want to implement at least 2 more vehicles, 2 more controllers, and 2 or 3 more estimators. And, finally, as far as the simulation is designed, I would like for the user to be able to select from a list of choices which vehicle, which controller, and which estimator to use. Therefore, I was looking for a clean way to keep a kind of registry of classes that implement each of these interfaces so that I wouldn't have to go back and change the interface code later when I implemented more sub-classes.

The solution that was accepted was to implement a registry class that maintains a mapping of class-names to constructors, and then update that mapping from within the class definition file from each of the implementing sub-classes using a static initializer. I went one step further and made the registry class generic (templated), and this is the result. There is an example code below.

File: CInterfaceRegistry.h

    /**
     *  file        CInterfaceRegistry.h
     *  date:      Aug 11, 2009
     *  brief:
     *
     *  detail:
     */

    #ifndef CINTERFACEREGISTRY_H_
    #define CINTERFACEREGISTRY_H_

    #include 
    #include 
    #include 
    #include 

    #include "exception/IllegalArgumentException.h"

    namespace utility
    {


    /**
     *  brief  Generic singleton object to maintains a registry of subclasses
     *          that implement a particular interface
     *
     *  This clever solution was taken from the discussion at the following site:
     *  http://stackoverflow.com/questions/1260954/how-can-i-enumerate-all-
     *  classes-that-implement-an-interface
     *
     *  note:   The registry will allow the registration of any class that extends
     *          an interface so long as it has a factory method, but the
     *          RegisterWithInterface macro will only work with classes that have
     *          named that method create()
     *
     *
     */
    template 
    class CInterfaceRegistry
    {
        private:
            std::map< std::string, InterfaceType* (*)(void) > m_creatorMap;


            /**
             *  brief  private to ensure that only the singleton ( the single
             *          static instantiation ) of this template class is used
             */
            CInterfaceRegistry(){}

        public:

            /**
             *  brief  returns the registry instance particular to the specified
             *          interface (specified as template parameter)
             */
            static CInterfaceRegistry& getInstance();


            /**
             *  brief  registers a new subclass of a an interface
             */
            bool registerClass( const std::string& name,
                                    InterfaceType* (*creator)(void) );


            /**
             *  brief  registers a new subclass from it's typeid() result, rather
             *          than from a hand-typed class name
             */
            bool registerClass( const std::type_info& classType,
                                    InterfaceType* (*creator)(void) );


            /**
             *  brief  returns a list of classes registered with this interface
             */
            std::set  getClassNames();


            /**
             *  brief  returns a new object of the specified class
             */
            InterfaceType*  createObjectOf( std::string className );

    };





    // A convient macro to compact the registration of a class
    #define RegisterWithInterface( CLASS, INTERFACE )                   
    namespace                                                           
    {                                                                   
        bool dummy_ ## CLASS =                                          
        CInterfaceRegistry::getInstance().registerClass(     
                typeid(CLASS), CLASS::create );                         
    }


    /**
     *  The use of this method ensures that this template class remains singleton.
     *  Only the static registry object created by this method will exist for any
     *  instantiation of this class.
     */
    template 
    CInterfaceRegistry& CInterfaceRegistry::getInstance()
    {
        static CInterfaceRegistry    registry;
        return registry;
    }


    /**
     *  To register a class with this registry, we map a factory function that is
     *  capable of generating objects of the class and returning a pointer of the
     *  interface type. The key we map it to is the name of the class.
     */
    template 
    bool CInterfaceRegistry::registerClass( const std::string& name,
                            InterfaceType* (*creator)(void) )
    {
        m_creatorMap[name] = creator;
        return true;
    }


    /**
     *  For added convenience, we can avoid typing the class names by hand if we
     *  use the typeid() operator. This method will extract the class name from
     *  a type_info class and pass it on to the actual registration method.
     */
    template 
    bool CInterfaceRegistry::registerClass( 
                            const std::type_info& classType,
                            InterfaceType* (*creator)(void) )
    {
        return registerClass( std::string(classType.name()), creator );
    }


    /**
     *  To generate a list of the class names currently registered with this
     *  interface, we iterate through the map and extract all the keys.
     */
    template 
    std::set  CInterfaceRegistry::getClassNames()
    {
        std::set   keys;

        typename
        std::map< std::string, InterfaceType* (*)(void) >::iterator pair;

        for( pair = m_creatorMap.begin(); pair != m_creatorMap.end(); pair++)
            keys.insert( pair->first );

        return keys;
    }


    /**
     *  To create a new object of the specified class, we simply de-reference the
     *  stored factory method and execute it.
     */
    template 
    InterfaceType*  CInterfaceRegistry::createObjectOf(
            std::string className )
    {
        InterfaceType* (*creator)(void) = m_creatorMap[className];

        if(creator)
            return *creator();
        else
            throw IllegalArgumentException(
                    className + "is not registered with the registry");
    }


    }

    #endif /* CINTERFACEREGISTRY_H_ */

File: CInterfaceRegistryTest.cpp

    /**
     *  file        CInterfaceRegistryTest.cpp
     *  date:      Aug 11, 2009
     *  brief:
     *
     *  detail:
     */

    #include "utility/CInterfaceRegistry.h"

    #include 
    #include 
    #include 

    using namespace utility;

    using std::map;
    using std::set;
    using std::cout;
    using std::string;



    // create the interface that we will be extending
    class IDummyInterface{};


    // create the first of several implementations of that interface
    class CDerivedA : public IDummyInterface
    {
        public:
            // the class must have some kind of static factory method
            static IDummyInterface* create(){ return new CDerivedA(); }
    };

    // and this is how we register the class with the registry
    bool dummyA =
            CInterfaceRegistry::getInstance().registerClass(
                    "CDerivedA", CDerivedA::create );


    // we create, here, the second of several implementations, it's basically the
    // same as the first
    class CDerivedB : public IDummyInterface
    {
        public:
            // again with the static factory method
            static IDummyInterface* create(){ return new CDerivedB(); }
    };

    // this is the same as above
    bool dummyB =
            CInterfaceRegistry::getInstance().registerClass(
                    "CDerivedB", CDerivedB::create );


    // and a nother implementation
    class CDerivedC : public IDummyInterface
    {
        public:
            // ditto…
            static IDummyInterface* create(){ return new CDerivedC(); }
    };

    // this time we use that convenient macro that does things a little more
    // compactly and, I think, without sacrificing readabilty
    RegisterWithInterface( CDerivedC, IDummyInterface );



    int main()
    {
        // here we can retrieve a list of all the registered classes by
        // querying the registry object
        set classes =
                CInterfaceRegistry::getInstance().getClassNames();

        cout << "Currently registered subclasses of IDummyInterface: n";
        cout << "----------------------nn";

        set::iterator   str;
        for( str = classes.begin(); str != classes.end(); str++ )
            cout << *str << "n";

        cout << "nndonen";

        return 0;
    }

Screen Capture of the test program

Note that the name returned by querying typeid() has a "9″ inserted in the front of it. The names used to identify a type in the type_info class are implementation specific, so it may or may not be a good idea to use them. In my case it will be fine.

Edit:
A better choice for the macro is to use something like this

    // A convient macro to compact the registration of a class
    #define RegisterWithInterface( CLASS, INTERFACE )                   
    namespace                                                           
    {                                                                   
        bool dummy_ ## CLASS =                                          
        Registry::getInstance().registerClass(     
                #CLASS, CLASS::createNew );                               
    }

which uses the class's name from the source code instead of the type_info name.

In-Place Matrix Transpose


I was looking around today for an algorithm for In-Place Matrix Transposition and didn't manage to find anything, so this is what I came up with. I'll analyze and document it fully later, but I just wanted to get it down.

It's wasteful in that it will repeat a lot of operations, especially for "thin" matrices (worst case: vector). But it will work, never-the-less. Worst case run time is at most $latex O(n^2)$ where $latex n$ is the size of the matrix $latex n = rows times columns$

Edit: No longer arbitrarily repeating cycles ( was leading to identity operation, duh! ).

    template 
    void CMatrix::inPlaceTranspose(  )
    {
        using std::swap;

        resize(cols,rows);

        int i, j, k, k_start, k_new;
        int n = rows;
        int m = cols;

        int length = rows*cols;

        for(k_start=1; k_start < length; k_start++)
        {
            T       temp    = data[k_start];
            bool    abort   = false;
            k_new = k = k_start;

            // go through the cycle once and ensure that the starting point is
            // minimal
            do
            {
                if( k_new < k_start )
                {
                    abort = true;
                    break;
                }

                k       = k_new;
                i       = k%n;
                j       = k/n;
                k_new   = i*m + j;
            }while(k_new != k_start);


            // if the current value is not the minimum of the cycle, then don't
            // perform the cycle
            if(abort)
                continue;


            // otherwise, perform the cycle
            k_new = k = k_start;
            do
            {
                swap( data[k_new], temp );

                k       = k_new;
                i       = k%n;
                j       = k/n;
                k_new   = i*m + j;
            }while(k_new != k_start);

            swap( data[k_new], temp );
        }
    }

And the result:
In Place Transpose