The points themselves belong to two subsets: for the first subset, let's call them "known" points, I know the coordinates from both images

Use regression to fit the parameters of the function transforming your points from the first subset to the second subset. Depending on the form of the distorsion, it could be multivariate linear regression, nonlinear regression you would need to provide a function (and starting values) for, or some sort of data-driven approach (P-splines, PLS, support vector machines with nonlinear basis functions, XGBoost...) that doesn't need an explicit function form and uses its many degrees of freedom to fit arbitrary function shapes. The PDL module family probably contains enough Levenberg-Marquardt and linear algebra to get you started.

and for the second subset, "target" points, I know their coordinates only from the reference image.

Use the parameters obtained above to perform the transformation.

If you didn't know the corresponding points in both images, the inner part of the problem would stay the same, and the outer part of the problem would use Iterative closest point to try to estimate and update those correspondences between points.