From 6c3a54027e9e8edde1efaa534371b7bac4519cec Mon Sep 17 00:00:00 2001
From: Saad Jbabdi <saad@fmrib.ox.ac.uk>
Date: Tue, 24 Nov 2009 18:37:42 +0000
Subject: [PATCH] added rotmat and rotwarp to deal with highres warpfields
---
vecreg.cc | 175 +++++++++++++++++++++++++++++++++++-------------------
1 file changed, 113 insertions(+), 62 deletions(-)
diff --git a/vecreg.cc b/vecreg.cc
index 0ea16f7..4b4f8d5 100755
--- a/vecreg.cc
+++ b/vecreg.cc
@@ -13,7 +13,7 @@
using namespace Utilities;
-string title="vecreg (Version 1.0)\nVector Affine/NonLinear Tranformation with Orientation Preservation";
+string title="vecreg (Version 1.1)\nVector Affine/NonLinear Tranformation with Orientation Preservation";
string examples="vecreg -i <input4D> -o <output4D> -r <refvol> -t <transform>";
Option<bool> verbose(string("-v,--verbose"),false,
@@ -37,8 +37,11 @@ Option<string> matrix(string("-t,--affine"),string(""),
Option<string> warp(string("-w,--warpfield"),string(""),
string("filename for 4D warp field for nonlinear registration"),
false,requires_argument);
-Option<string> matrix2(string("--postmat"),string(""),
- string("filename for post-transform affine matrix \n\t\t\tif set, this matrix will be used for the rotation of the vector/tensor field "),
+Option<string> matrix2(string("--rotmat"),string(""),
+ string("filename for secondary affine matrix \n\t\t\tif set, this will be used for the rotation of the vector/tensor field "),
+ false,requires_argument);
+Option<string> warp2(string("--rotwarp"),string(""),
+ string("filename for secondary warp field \n\t\t\tif set, this will be used for the rotation of the vector/tensor field "),
false,requires_argument);
Option<string> interpmethod(string("--interp"),"",
string("interpolation method : nearestneighbour, trilinear (default), sinc or spline"),
@@ -46,6 +49,9 @@ Option<string> interpmethod(string("--interp"),"",
Option<string> maskfile(string("-m,--mask"),string(""),
string("brain mask in input space"),
false,requires_argument);
+Option<string> omaskfile(string("--refmask"),string(""),
+ string("brain mask in output space (useful for speed up of nonlinear reg)"),
+ false,requires_argument);
////////////////////////////////////////////////////////
ReturnMatrix rodrigues(const float& angle,ColumnVector& w){
@@ -137,6 +143,33 @@ ReturnMatrix ppd(const Matrix& F,ColumnVector& e1){
return R;
}
+void sjgradient(const volume<float>& im,volume4D<float>& grad){
+
+ grad.reinitialize(im.xsize(),im.ysize(),im.zsize(),3);
+ copybasicproperties(im,grad[0]);
+
+ int fx,fy,fz,bx,by,bz;
+ float dx,dy,dz;
+ for (int z=0; z<grad.zsize(); z++){
+ fz = z ==(grad.zsize()-1) ? 0 : 1;
+ bz = z == 0 ? 0 : -1;
+ dz = (fz==0 || bz==0) ? 1.0 : 2.0;
+ for (int y=0; y<grad.ysize(); y++){
+ fy = y ==(grad.ysize()-1) ? 0 : 1;
+ by = y == 0 ? 0 : -1;
+ dy = (fy==0 || by==0) ? 1.0 : 2.0;
+ for (int x=0; x<grad.xsize(); x++){
+ fx = x ==(grad.xsize()-1) ? 0 : 1;
+ bx = x == 0 ? 0 : -1;
+ dx = (fx==0 || bx==0) ? 1.0 : 2.0;
+ grad[0](x,y,z) = (im(x+fx,y,z) - im(x+bx,y,z))/dx;
+ grad[1](x,y,z) = (im(x,y+fy,z) - im(x,y+by,z))/dy;
+ grad[2](x,y,z) = (im(x,y,z+fz) - im(x,y,z+bz))/dz;
+ }
+ }
+ }
+
+}
void vecreg_aff(const volume4D<float>& tens,
volume4D<float>& oV1,
@@ -147,12 +180,26 @@ void vecreg_aff(const volume4D<float>& tens,
Matrix iM(4,4);
iM=M.i();
- //////////////////////////////////////////////////
- // Do we only need a FLIRT matrix for reorienting?
+ /////////////////////////////////////////////////////////////////////////
+ // Where we define a potential affine transfo for the rotation
Matrix M2;
if(matrix2.value()!="")
M2 = read_ascii_matrix(matrix2.value());
- //////////////////////////////////////////////////
+ // Where we define a potential warp field transfo for the rotation
+ volume4D<float> warpvol2;
+ volume4D<float> jx,jy,jz;
+ Matrix Jw(3,3),I(3,3);I<<1<<0<<0<<0<<1<<0<<0<<0<<1;
+ if(warp2.value()!=""){
+ FnirtFileReader ffr2(warp2.value());
+ warpvol2=ffr2.FieldAsNewimageVolume4D(true);
+ sjgradient(warpvol2[0],jx);
+ sjgradient(warpvol2[1],jy);
+ sjgradient(warpvol2[2],jz);
+ }
+ /////////////////////////////////////////////////////////////////////////
+ volume<float> omask;
+ if(omaskfile.value()!="")
+ read_volume(omask,omaskfile.value());
// extract rotation matrix from M
@@ -174,14 +221,17 @@ void vecreg_aff(const volume4D<float>& tens,
for(int z=0;z<oV1.zsize();z++)
for(int y=0;y<oV1.ysize();y++)
for(int x=0;x<oV1.xsize();x++){
+ if(omaskfile.value()!="")
+ if(omask(x,y,z)==0)
+ continue;
// compute seed coordinates
X_target << x << y << z;
X_seed=vox_to_vox(X_target,targetdim,seeddim,iM);
- if(mask((int)X_seed(1),(int)X_seed(2),(int)X_seed(3))==0){
+ if(mask((int)round(float(X_seed(1))),(int)round(float(X_seed(2))),(int)round(float(X_seed(3))))==0)
continue;
- }
+
// compute interpolated tensor
Tens.Row(1) << tens[0].interpolate(X_seed(1),X_seed(2),X_seed(3));
@@ -192,11 +242,21 @@ void vecreg_aff(const volume4D<float>& tens,
<< tens[5].interpolate(X_seed(1),X_seed(2),X_seed(3));
+
+ if(warp2.value()!=""){
+ // Local Jacobian of the backward warpfield
+ Jw << jx(x,y,z,0) << jx(x,y,z,1) << jx(x,y,z,2)
+ << jy(x,y,z,0) << jy(x,y,z,1) << jy(x,y,z,2)
+ << jz(x,y,z,0) << jz(x,y,z,1) << jz(x,y,z,2);
+ // compute local forward affine transformation
+ F = (I + Jw).i();
+ }
+
if(oV1.tsize()==3){ // case where input is a vector
// compute first eigenvector
EigenValues(Tens,d,v);
V_seed = v.Column(3);
-
+
// rotate vector
V_target=F*V_seed;
if(V_target.MaximumAbsoluteValue()>0)
@@ -209,7 +269,14 @@ void vecreg_aff(const volume4D<float>& tens,
// create tensor
if(oV1.tsize()==6){
- Tens << R*Tens*R.t();
+ if(warp2.value()!=""){
+ EigenValues(Tens,d,v);
+ R=ppd(F,v.Column(3),v.Column(2));
+ Tens << R*Tens*R.t();
+ }
+ else{
+ Tens << R*Tens*R.t();
+ }
oV1(x,y,z,0)=Tens(1,1);
oV1(x,y,z,1)=Tens(2,1);
@@ -223,33 +290,6 @@ void vecreg_aff(const volume4D<float>& tens,
}
}
-void sjgradient(const volume<float>& im,volume4D<float>& grad){
-
- grad.reinitialize(im.xsize(),im.ysize(),im.zsize(),3);
- copybasicproperties(im,grad[0]);
-
- int fx,fy,fz,bx,by,bz;
- float dx,dy,dz;
- for (int z=0; z<grad.zsize(); z++){
- fz = z ==(grad.zsize()-1) ? 0 : 1;
- bz = z == 0 ? 0 : -1;
- dz = (fz==0 || bz==0) ? 1.0 : 2.0;
- for (int y=0; y<grad.ysize(); y++){
- fy = y ==(grad.ysize()-1) ? 0 : 1;
- by = y == 0 ? 0 : -1;
- dy = (fy==0 || by==0) ? 1.0 : 2.0;
- for (int x=0; x<grad.xsize(); x++){
- fx = x ==(grad.xsize()-1) ? 0 : 1;
- bx = x == 0 ? 0 : -1;
- dx = (fx==0 || bx==0) ? 1.0 : 2.0;
- grad[0](x,y,z) = (im(x+fx,y,z) - im(x+bx,y,z))/dx;
- grad[1](x,y,z) = (im(x,y+fy,z) - im(x,y+by,z))/dy;
- grad[2](x,y,z) = (im(x,y,z+fz) - im(x,y,z+bz))/dz;
- }
- }
- }
-
-}
void vecreg_nonlin(const volume4D<float>& tens,volume4D<float>& oV1,
@@ -258,11 +298,6 @@ void vecreg_nonlin(const volume4D<float>& tens,volume4D<float>& oV1,
ColumnVector X_seed(3),X_target(3);
- //float dxx=tens.xdim(),dyy=tens.ydim(),dzz=tens.zdim();
- //float dx=oV1.xdim(),dy=oV1.ydim(),dz=oV1.zdim();
- //float nxx=(float)tens.xsize()/2.0,nyy=(float)tens.ysize()/2.0,nzz=(float)tens.zsize()/2.0;
- //float nx=(float)oV1.xsize()/2.0,ny=(float)oV1.ysize()/2.0,nz=(float)oV1.zsize()/2.0;
-
// read warp field created by Jesper
FnirtFileReader ffr(warp.value());
@@ -271,8 +306,8 @@ void vecreg_nonlin(const volume4D<float>& tens,volume4D<float>& oV1,
Matrix F(3,3),u(3,3),v(3,3);
DiagonalMatrix d(3);
- //////////////////////////////////////////////////
- // Do we only need a FLIRT matrix for reorienting?
+ /////////////////////////////////////////////////////////////////////////
+ // Where we define a potential affine transfo for the rotation
Matrix M2;
if(matrix2.value()!=""){
M2 = read_ascii_matrix(matrix2.value());
@@ -281,16 +316,25 @@ void vecreg_nonlin(const volume4D<float>& tens,volume4D<float>& oV1,
SVD(F*F.t(),d,u,v);
F=(u*sqrt(d)*v.t()).i()*F;
}
- //////////////////////////////////////////////////
-
- // compute transformation jacobian
- volume4D<float> jx(mask.xsize(),mask.ysize(),mask.zsize(),3);
- volume4D<float> jy(mask.xsize(),mask.ysize(),mask.zsize(),3);
- volume4D<float> jz(mask.xsize(),mask.ysize(),mask.zsize(),3);
- sjgradient(warpvol[0],jx);
- sjgradient(warpvol[1],jy);
- sjgradient(warpvol[2],jz);
-
+ // Where we define a potential warp field transfo for the rotation
+ volume4D<float> warpvol2;
+ volume4D<float> jx,jy,jz;
+ if(warp2.value()!=""){
+ FnirtFileReader ffr2(warp2.value());
+ warpvol2=ffr2.FieldAsNewimageVolume4D(true);
+ sjgradient(warpvol2[0],jx);
+ sjgradient(warpvol2[1],jy);
+ sjgradient(warpvol2[2],jz);
+ }
+ else{
+ sjgradient(warpvol[0],jx);
+ sjgradient(warpvol[1],jy);
+ sjgradient(warpvol[2],jz);
+ }
+ /////////////////////////////////////////////////////////////////////////
+ volume<float> omask;
+ if(omaskfile.value()!="")
+ read_volume(omask,omaskfile.value());
ColumnVector V_seed(3),V_target(3);
ColumnVector V1_seed(3),V2_seed(3);
@@ -301,16 +345,16 @@ void vecreg_nonlin(const volume4D<float>& tens,volume4D<float>& oV1,
for(int z=0;z<oV1.zsize();z++)
for(int y=0;y<oV1.ysize();y++)
for(int x=0;x<oV1.xsize();x++){
-
-
+ if(omaskfile.value()!="")
+ if(omask(x,y,z)==0)
+ continue;
X_target << x << y << z;
X_seed = NewimageCoord2NewimageCoord(warpvol,false,oV1[0],mask,X_target);
- if(mask((int)X_seed(1),(int)X_seed(2),(int)X_seed(3))==0){
- continue;
- }
+ if(mask((int)round(float(X_seed(1))),(int)round(float(X_seed(2))),(int)round(float(X_seed(3))))==0)
+ continue;
// compute interpolated tensor
Tens.Row(1) << tens[0].interpolate(X_seed(1),X_seed(2),X_seed(3));
@@ -448,6 +492,7 @@ int do_vecreg(){
}
//cout<<"elapsed time:"<<time(NULL)-_time<<" sec"<<endl;
+ ovol.setDisplayMaximumMinimum(ivol.max(),ivol.min());
save_volume4D(ovol,outfilename.value());
return 0;
@@ -467,10 +512,12 @@ int main(int argc,char *argv[]){
options.add(outfilename);
options.add(ref);
options.add(matrix);
- options.add(matrix2);
options.add(warp);
+ options.add(matrix2);
+ options.add(warp2);
options.add(interpmethod);
options.add(maskfile);
+ options.add(omaskfile);
options.parse_command_line(argc,argv);
@@ -480,19 +527,23 @@ int main(int argc,char *argv[]){
exit(EXIT_FAILURE);
}
if( (matrix.set()) && (warp.set()) ){
- options.usage();
cerr << endl
<< "Cannot specify both --affine AND --warpfield"
<< endl << endl;
exit(EXIT_FAILURE);
}
if( (matrix.unset()) && (warp.unset()) ){
- options.usage();
cerr << endl
<< "Please Specify either --affine OR --warpfield"
<< endl << endl;
exit(EXIT_FAILURE);
}
+ if( (warp2.set()) && (matrix2.set()) ){
+ cerr << endl
+ << "Cannot Specify both --rotaff AND --rotwarp"
+ << endl << endl;
+ exit(EXIT_FAILURE);
+ }
}
catch(X_OptionError& e) {
--
GitLab