Make a quad, position at midpoint of p1 and p2. (trivial)
Scale it on lets say X by distance between p1 and p2. So the position matrix looks like:
(p1-p2).magnitude(), 0 , 0 , 0
0 , 1 , 0 , 0
0 , 0 , 1 , 0
midpoint.x , midpoint.y, midpoint.z, 1
(thats row major i think. bottom row is position)
Now you just need the matrix that looks from billboard to camera and use some math to combine them (multiply? in my own matrix code i need to do rotation and translation separately for it to work but it was written by me so... yeah.)
For the rotation matrix youll need:
-X axis direction (direction from midpoint/p1 to p2)
-Y axis direction (cross the above and the below vectors (X axis dir and Z axis dir))
-Z axis direction (direction from midpoint to camera)
So now we only need the direction from midpoint to camera, which is simply taking the unit vector from their difference.
What you need to know:
[local X axis direction] ,0
[local Y axis direction] ,0
[local Z axis direction] ,0
[position ] , 1 //1 means position 0 means vector/direction
Takes 2 vectors, gives a vector that is perpendicular to both (imagine the 2 vectors form a plane. The result is the normal of the plane. The 2 input vectors dont need to be perpendicular with each other!)
sqrt(xDifference^2 + yDifference^2 + zDifference^2)
Take the difference of 2 vectors (a-b), divide by the magnitude of the result. Gives direction vector with length 1.
That should be all that you need.