Just a few words about a work I am doing for G’MIC.
I’ve started to implement the Iterative Closest Point algorithm in G’MIC. This algorithm can be used to register two point clouds together (so basically two shapes as well) :
What I’ve done so far is just testing stuffs, so this is only the beginning. But my first tests are promising:
- First experiment : Register a star shape (in red) and a rotated/de-zoomed star shape (in green).
(red shape is the target, green shape is the source and blue lines correspond to the closest matchs found at each iteration). The source / target shapes do not have the same number of points, of course.
- Second experiment : Register a heart shape with a noisy heart shape. Even with the fact the source shape (green) has noisy coordinates, the algorithm converges into something correct.
The transformation here has 6 parameters and is a mix of rotation/translation/shearing.
As I said, this is only a prototype for now, but I’d like to make it available as a stdlib command in the future. I’m pretty sure we could find interesting uses of this
(@patdavid , you know what I’m talking about! 
).
Feel free to play with the current prototype code if you wish :
foo :
shape_heart 300 +erode. 3 -[-2,-1]
1,1,1,2 eval.. "i?da_push([x,y]); end(resize(#-1,1,da_size(),1,2,0))" rm..
shape_heart 200 rotate. 35 b. 1 ge. 50% +erode. 3 -[-2,-1] r. 300,300,1,1,0,0,0.5,0.5
1,1,1,2 eval.. "i?da_push([x,y]); end(resize(#-1,1,da_size(),1,2,0))" rm..
# Start iteration loop.
repeat inf,f
# Find closest points
1,1,1,2 eval.. ">
X = [i0,i1];
dmin = inf;
Pmin = [ 0,0 ];
repeat (h#0,p,
P = I[#0,p];
d = norm(P - X);
d<dmin?(dmin = d; Pmin = P);
);
critical(da_push(Pmin));
I;
end(resize(#-1,1,da_size(),1,2,0));
"
300,300,1,3
eval[0] ">I(#-1,i0,i1) = [ 255,0,0 ]"
eval[1] ">P = [ i0,i1]; Q = I[#2,y]; polygon(#-1,2,P,Q,0.5,0,0,255)"
eval[1] ">I(#-1,i0,i1) = [ 0,255,0 ]"
w. 600,600
rm.
# Solve least-square problem.
+permute[2] cyzx y. # Y
+permute[1] cyzx 1,100%,1,1,1 a[-2,-1] x
l. r. 100%,200%,1,1,1 s y repeat $! r[$>] 6,100%,1,1,0,0,{$>%2} done a y endl # A
solve.. . rm.
rm..
# Apply estimated transformation on source coordinates.
f[1] "
begin(
A = i[#-1,0];
B = i[#-1,1];
C = i[#-1,2];
D = i[#-1,3];
E = i[#-1,4];
F = i[#-1,5];
);
X = i0;
Y = i1;
[ A*X + B*Y + C,
D*X + E*Y + F ]"
rm.
done
