Mat4x4     skew;
     Mat4x4 dskew[3];
     
     Vec3d axisUn(rodParms[0], rodParms[1], rodParms[2]);
     Vec3d axis = axisUn;
     Vec3d dangle(axis.x, axis.y, axis.z);
         
     double angle = std::max(axis.normalize(), 1e-5);
     double sa = sin(angle), ca = cos(angle);
 
     dangle *= (1.0/angle);
 
     memset(skew.data, 0, sizeof(Mat4x4));
     memset(dskew, 0, sizeof(Mat4x4)*3);
 
     // Axis is normalized.
     skew(1, 2) = -axis.x;
     skew(2, 1) =  axis.x;
 
     skew(2, 0) = -axis.y;
     skew(0, 2) =  axis.y;
 
     skew(0, 1) = -axis.z;
     skew(1, 0) =  axis.z;
        
     /*
       (d/dx)rx/angle ==> (angle * ddx rx - rx * ddx angle)/(angle*angle)
     */
     for(int k=0; k3; k++){
       dskew[k](1, 2) = -(angle * (k==0) - rodParms[0] * dangle.data[k])/(angle*angle);
       dskew[k](2, 1) =  (angle * (k==0) - rodParms[0] * dangle.data[k])/(angle*angle);
 
       dskew[k](2, 0) = -(angle * (k==1) - rodParms[1] * dangle.data[k])/(angle*angle);
       dskew[k](0, 2) =  (angle * (k==1) - rodParms[1] * dangle.data[k])/(angle*angle);
 
       dskew[k](0, 1) = -(angle * (k==2) - rodParms[2] * dangle.data[k])/(angle*angle);
       dskew[k](1, 0) =  (angle * (k==2) - rodParms[2] * dangle.data[k])/(angle*angle);
     }
 
     Mat4x4 trans = Mat4x4::trans(rodParms[3], rodParms[4], rodParms[5]);
     Mat4x4 rot = Mat4x4::eye() + skew * sa + (skew * skew)*(1 - ca);
 
     memset(derivs, 0, sizeof(Mat4x4)*6);
 
     for(int k=0; k3; k++)
       derivs[k] = trans*(dskew[k]*sa + skew*(ca*dangle[k]) 
                        + (skew*dskew[k] + dskew[k]*skew)*(1-ca) + skew*skew*(sa*dangle[k]));                
     derivs[3](0, 3) = 1;
     derivs[4](1, 3) = 1;
     derivs[5](2, 3) = 1;