simulate.m

#!/bin/octave
## @ moritz siegel

clear all
close all
clc
graphics_toolkit('gnuplot')
#graphics_toolkit('qt')
#graphics_toolkit('fltk')

## \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
## soft settings \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
global hd = "~/biophysics"; # parent directory
global wd = "data";
global nn = "AFFINE"; # additional luminosity suffix  ( default: "" )
method = "affine"; # /char, "rigid","affine", method to reconstruct
#method = "affine"; # /char, "rigid","affine", method to reconstruct
eps = 1e-8; # precision ( default: 1e-8 )
## end of settings: hands off /////////////////////////////////////////////////
## ///////////////////////////////////////////////////////////////////////////


disp( 'init \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ ' );

## lets move it move it.
addpath( hd ) # function & stuff needed
stamp = strftime("%Y_%m_%d_%H%M%S", localtime (time ())); # create timestamp for saving files
global nwd = sprintf( "%s/%s/simulation_%s_%s_%s", hd, wd, method, stamp, nn ); # concat working dir
mkdir( nwd )
chdir( nwd )

## init first point cloud. introduce noise for z component, matrix cant
## be singular (e.g. points on a plane), else cholesky decomposition fails.
n = 100;
theta = 100 * rand( n, 1 );
p = [ cos( theta ), sin( theta ), 1e-3 * rand( size( theta ) ) ];

## introduce noise
salt = [ 1e-2*rand( size( theta ) ), 1e-2*rand( size( theta ) ), 1e-5*rand( size( theta ) ) ];
q = p + salt;

## define affine transforms & derive second point cloud.
simulated_translation = [ 0, 0.005, 0.003 ];
reflect = [ -1, 0, 0; 0, 1, 0; 0, 0, 1 ];
scale = [ 2, 0, 0; 0, 1, 0; 0, 0, 1 ];
theta = 70;
rot = [ cosd(theta), 0, -sind(theta); 0, 1, 0; sind(theta), 0, cosd(theta) ];
shear = [ 1, 0.5, 0; 0, 1, 0; 0, 0, 1 ];
simulated_transform = eye( 3 );
simulated_transform = simulated_transform * rot;
simulated_transform = simulated_transform * reflect;
simulated_transform = simulated_transform * scale;
simulated_transform = simulated_transform * shear;
q = ( simulated_transform * q' )';
q = q + simulated_translation;

## check determinant.
#assert( det( simulated_transform ) > 0, "det( simulated_transform ) <= 0, rerun" );


disp( 'analyse \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ ' );

## recover transform
switch( method )
case "rigid"
    [ rekovered_transform, rekovered_translation ] = rig( p, q );
case "affine"
    [ rekovered_transform, rekovered_translation ] = affine( p, q );
endswitch


disp( 'simulated translation vector:' )
disp( simulated_translation )
disp( 'recovered translation vector:' )
disp( rekovered_translation )
disp( 'deviation translation vector:' )
disp( simulated_translation - rekovered_translation )
disp( 'simulated affine transform:' )
disp( simulated_transform )
disp( 'recovered affine transform:' )
disp( rekovered_transform )
disp( 'deviation affine transform:' )
disp( simulated_transform - rekovered_transform )

save simulated_translation.mat simulated_translation
save rekovered_translation.mat rekovered_translation
save simulated_transform.mat simulated_transform
save rekovered_transform.mat rekovered_transform

if ( any( any( ( simulated_transform - rekovered_transform ) > eps ) ) )
  disp( "warning: failed to recover affine transform" );
endif


disp( 'plot \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ ' );

## compute translation.
p_recovered  = ( rekovered_transform * p' )' + rekovered_translation;

## plot & compare both orthonormalised point clouds with shifted.
fhor0 = figure;
plot3( p(:,1), p(:,2), p(:,3), '.b' );
hold on;
axis equal;
plot3( q(:,1), q(:,2), q(:,3), '.r' );
plot3( p_recovered(:,1), p_recovered(:,2), p_recovered(:,3), 'b' );
print( fhor0, "affine_5_affine_rotated_vs_shifted.png" );

dlmwrite( "netz.dat", "x y z", " " )
dlmwrite( "p.dat", "x y z", " " )
dlmwrite( "q.dat", "x y z", " " )
dlmwrite( "pr.dat", "x y z", " " )
for c = 1:length( p )
    dlmwrite( "p.dat", p(c,:), " ", "-append" )
    dlmwrite( "q.dat", q(c,:), " ", "-append" )
    dlmwrite( "pr.dat", p_recovered(c,:), " ", "-append" )
    dlmwrite( "netz.dat", p(c,:), " ", "-append" )
    dlmwrite( "netz.dat", p_recovered(c,:), " ", "-append" )
    dlmwrite( "netz.dat", "", " ", "-append" )
endfor