some doubts about arGetTransMatSub() function of artoolkit in the 3d registration step

Question

I'm learning artoolkit codes, but encountered some difficulties, it has cost me more than two weeks, especially in the three dimensional registration part, I don't really know arGetTransMatSub () this function, especially the inside of the algorithm and variable. I've read some documents about artoolkit three dimensional registration , but I can't put those theory and codes combined, the following are some of my questions:
1. Is the function of ArParamIdeal2Observ () is to translate the coordinates of marker in the Ideal Screen Coordinates system to the Observed Screen Coordinates system?
2. In the for loop block, is wx, wy, wz represents the coordinates of the marker in the camera coordinate system?
3.What mat_b, mat_c, mat_d , mat_e respectively represents for?
I really want to master these, I'm appreciated for your help.

Below is some codes from arGetTransMatSub() function.

if( arFittingMode == AR_FITTING_TO_INPUT ) 
{
    for( i = 0; i < num; i++ ) 
    {
        arParamIdeal2Observ(dist_factor, ppos2d[i][0], ppos2d[i][1],
                                         &pos2d[i][0], &pos2d[i][1]);
    }
}
else
{
    for( i = 0; i < num; i++ ) 
    {
        pos2d[i][0] = ppos2d[i][0];
        pos2d[i][1] = ppos2d[i][1];
    }
}

for( j = 0; j < num; j++ ) 
{
    wx = rot[0][0] * pos3d[j][0]
       + rot[0][1] * pos3d[j][1]
       + rot[0][2] * pos3d[j][2];
    wy = rot[1][0] * pos3d[j][0]
       + rot[1][1] * pos3d[j][1]
       + rot[1][2] * pos3d[j][2];
    wz = rot[2][0] * pos3d[j][0]
       + rot[2][1] * pos3d[j][1]
       + rot[2][2] * pos3d[j][2];
    mat_a->m[j*6+0] = mat_b->m[num*0+j*2] = cpara[0][0];
    mat_a->m[j*6+1] = mat_b->m[num*2+j*2] = cpara[0][1];
    mat_a->m[j*6+2] = mat_b->m[num*4+j*2] = cpara[0][2] - pos2d[j][0];
    mat_c->m[j*2+0] = wz * pos2d[j][0]
           - cpara[0][0]*wx - cpara[0][1]*wy - cpara[0][2]*wz;
    mat_a->m[j*6+3] = mat_b->m[num*0+j*2+1] = 0.0;
    mat_a->m[j*6+4] = mat_b->m[num*2+j*2+1] = cpara[1][1];
    mat_a->m[j*6+5] = mat_b->m[num*4+j*2+1] = cpara[1][2] - pos2d[j][1];
    mat_c->m[j*2+1] = wz * pos2d[j][1]
          - cpara[1][1]*wy - cpara[1][2]*wz;
}
arMatrixMul( mat_d, mat_b, mat_a );
arMatrixMul( mat_e, mat_b, mat_c );
arMatrixSelfInv( mat_d );
arMatrixMul( mat_f, mat_d, mat_e );
trans[0] = mat_f->m[0];
trans[1] = mat_f->m[1];
trans[2] = mat_f->m[2];

ret = arModifyMatrix( rot, trans, cpara, pos3d, pos2d, num );