Calculation of 4x4 matrix of planes

Based on the code of vtkImageResliceMapper.cxx, I created a program to calculate 4 x 4 matrix on a plane.

The code created is shown below.

  // Setup a visualization pipeline
  vtkSmartPointer<vtkPlane> plane =vtkSmartPointer<vtkPlane>::New();
  double normal[3];
  plane->GetNormal(normal);
  double Origin[3];
  plane->GetOrigin(Origin);
  // Set the slice orientation
 vtkSmartPointer<vtkMatrix4x4> resliceMatrix =
vtkSmartPointer<vtkMatrix4x4>::New();

// Create the corresponding axis
double axis[3];
axis[0] = 0.0;
axis[1] = 0.0;
axis[2] = 0.0;

double saxis[3], taxis[3];
taxis[0] = 0.0;
taxis[1] = 1.0;
taxis[2] = 0.0;
vtkMath::Cross(taxis, axis, saxis);

// Compute the rotation angle between the axis and the normal
double vec[3];
vtkMath::Cross(axis, normal, vec);
double costheta = vtkMath::Dot(axis, normal);
double sintheta = vtkMath::Norm(vec);
double theta = atan2(sintheta, costheta);
if (sintheta != 0)
{
  vec[0] /= sintheta;
  vec[1] /= sintheta;
  vec[2] /= sintheta;
}

// convert to quaternion
costheta = cos(0.5*theta);
sintheta = sin(0.5*theta);
double quat[4];
quat[0] = costheta;
quat[1] = vec[0]*sintheta;
quat[2] = vec[1]*sintheta;
quat[3] = vec[2]*sintheta;

// convert to matrix
double mat[3][3];
vtkMath::QuaternionToMatrix3x3(quat, mat);

// Create a slice-to-data transform matrix
// The columns are v1, v2, normal
double v1[3], v2[3];
vtkMath::Multiply3x3(mat, saxis, v1);
vtkMath::Multiply3x3(mat, taxis, v2);

resliceMatrix->Element[0][0] = v1[0];
resliceMatrix->Element[1][0] = v1[1];
resliceMatrix->Element[2][0] = v1[2];
resliceMatrix->Element[3][0] = 0.0;

resliceMatrix->Element[0][1] = v2[0];
resliceMatrix->Element[1][1] = v2[1];
resliceMatrix->Element[2][1] = v2[2];
resliceMatrix->Element[3][1] = 0.0;

resliceMatrix->Element[0][2] = normal[0];
resliceMatrix->Element[1][2] = normal[1];
resliceMatrix->Element[2][2] = normal[2];
resliceMatrix->Element[3][2] = 0.0;

resliceMatrix->Element[0][3] = Origin[0];
resliceMatrix->Element[1][3] = Origin[1];
resliceMatrix->Element[2][3] = Origin[2];
resliceMatrix->Element[3][3] = 1.0;

Can you tell me if there is something to be improved?

@dgobbi