diff --git a/miscmaths.cc b/miscmaths.cc
index 295bb89a0386379c336186930408f02314aa90fc..aa7b78d9afac164bb06d42a0ca778873062851d1 100644
--- a/miscmaths.cc
+++ b/miscmaths.cc
@@ -560,9 +560,9 @@ namespace MISCMATHS {
   ReturnMatrix sqrtaff(const Matrix& mat)
     {
       Tracer tr("sqrtaff");
-      Matrix matnew(4,4), rot(4,4), id4(4,4);
-      Identity(rot);
-      Identity(id4);
+      Matrix matnew(4,4), rot, id4;
+      rot=IdentityMatrix(4);
+      id4=IdentityMatrix(4);
       ColumnVector params(12), centre(3), trans(4);
       centre = 0.0;
       // Quaternion decomposition -> params(1..3) = sin(theta/2)*(unit_axis_vec)
@@ -591,14 +591,12 @@ namespace MISCMATHS {
       rot(2,4) = 0.0;
       rot(3,4) = 0.0;
       
-      Matrix scale(4,4);
-      Identity(scale);
+      Matrix scale=IdentityMatrix(4);
       scale(1,1)=params(7);
       scale(2,2)=params(8);
       scale(3,3)=params(9);
 
-      Matrix skew(4,4);
-      Identity(skew);
+      Matrix skew=IdentityMatrix(4);
       skew(1,2)=params(10);
       skew(1,3)=params(11);
       skew(2,3)=params(12);
@@ -694,7 +692,7 @@ namespace MISCMATHS {
       Tracer tr("construct_rotmat_euler");
       ColumnVector angl(3);
       Matrix newaff(4,4);
-      Identity(aff);
+      aff=IdentityMatrix(4);
 
       if (n<=0) return 0;
       // order of parameters is 3 rotation + 3 translation
@@ -741,7 +739,7 @@ namespace MISCMATHS {
 			    const ColumnVector& centre)
     {
       Tracer tr("construct_rotmat_quat");
-      Identity(aff);
+      aff=IdentityMatrix(4);
 
       if (n<=0) return 0;
       // order of parameters is 3 rotation (last 3 quaternion components) 
@@ -793,7 +791,7 @@ namespace MISCMATHS {
     {
       // Matrix rot must be 4x4; angl and orig must be length 3
       Tracer tr("make_rot");
-      Identity(rot);  // default return value
+      rot=IdentityMatrix(4);  // default return value
       float theta;
       theta = norm2(angl);
       if (theta<1e-8) {  // avoid round-off errors and return Identity
@@ -815,8 +813,7 @@ namespace MISCMATHS {
       basischange.SubMatrix(1,3,2,2) = x3;
       basischange.SubMatrix(1,3,3,3) = x1;
 
-      Matrix rotcore(3,3);
-      Identity(rotcore);
+      Matrix rotcore=IdentityMatrix(3);
       rotcore(1,1)=cos(theta);
       rotcore(2,2)=cos(theta);
       rotcore(1,2)=sin(theta);
@@ -824,8 +821,7 @@ namespace MISCMATHS {
 
       rot.SubMatrix(1,3,1,3) = basischange * rotcore * basischange.t();
   
-      Matrix ident3(3,3);
-      Identity(ident3);
+      Matrix ident3=IdentityMatrix(3);
       ColumnVector trans(3);
       trans = (ident3 - rot.SubMatrix(1,3,1,3))*centre;
       rot.SubMatrix(1,3,4,4)=trans;
@@ -837,12 +833,12 @@ namespace MISCMATHS {
     {
       Tracer tr("getrotaxis");
       Matrix residuals(3,3);
-      residuals = rotmat*rotmat.t() - Identity(3);
+      residuals = rotmat*rotmat.t() - IdentityMatrix(3);
       if (residuals.SumSquare() > 1e-4)
 	{ cerr << "Failed orthogonality check!" << endl;  return -1; }
       Matrix u(3,3), v(3,3);
       DiagonalMatrix d(3);
-      SVD(rotmat-Identity(3),d,u,v);
+      SVD(rotmat-IdentityMatrix(3),d,u,v);
       // return column of V corresponding to minimum value of |S|
       for (int i=1; i<=3; i++) {
 	if (fabs(d(i))<1e-4)  axis = v.SubMatrix(1,3,i,i);
@@ -994,8 +990,7 @@ namespace MISCMATHS {
   
       if (n<=6)  return 0;
   
-      Matrix scale(4,4);
-      Identity(scale);
+      Matrix scale=IdentityMatrix(4);
       if (n>=7) {
 	scale(1,1)=params(7);
 	if (n>=8) scale(2,2)=params(8);
@@ -1008,8 +1003,7 @@ namespace MISCMATHS {
       strans = centre - scale.SubMatrix(1,3,1,3)*centre;
       scale.SubMatrix(1,3,4,4) = strans;
 
-      Matrix skew(4,4);
-      Identity(skew);
+      Matrix skew=IdentityMatrix(4);
       if (n>=10) {
 	if (n>=10) skew(1,2)=params(10);
 	if (n>=11) skew(1,3)=params(11);
@@ -1032,7 +1026,7 @@ float rms_deviation(const Matrix& affmat1, const Matrix& affmat2,
   Tracer trcr("rms_deviation");
   Matrix isodiff(4,4);
   try {
-    isodiff = affmat1*affmat2.i() - Identity(4);
+    isodiff = affmat1*affmat2.i() - IdentityMatrix(4);
   } catch(...) {
     cerr << "RMS_DEVIATION ERROR:: Could not invert matrix" << endl;  
     exit(-5); 
@@ -1091,7 +1085,7 @@ void get_axis_orientations(const Matrix& sform_mat, int sform_code,
     vox2mm = qform_mat;
   } else {
     // ideally should be sampling_mat(), but for orientation it doesn't matter
-    vox2mm = Identity(4);
+    vox2mm = IdentityMatrix(4);
     vox2mm(1,1) = -vox2mm(1,1);
   }
   mat44 v2mm;
@@ -1977,7 +1971,7 @@ void detrend(Matrix& p_ts, int p_level)
     }
           
   // Form residual forming matrix R:
-  Matrix R = Identity(sizeTS)-a*pinv(a);
+  Matrix R = IdentityMatrix(sizeTS)-a*pinv(a);
 
   for(int t = 1; t <= sizeTS; t++)
     {
@@ -2051,7 +2045,7 @@ void ols(const Matrix& data,const Matrix& des,const Matrix& tc, Matrix& cope,Mat
   }  
   Matrix pdes = pinv(des);
   Matrix prevar=diag(tc*pdes*pdes.t()*tc.t());
-  Matrix R=Identity(des.Nrows())-des*pdes;
+  Matrix R=IdentityMatrix(des.Nrows())-des*pdes;
   float tR=R.Trace();
   Matrix pe=pdes*data;
   cope=tc*pe;
@@ -2064,7 +2058,7 @@ void ols(const Matrix& data,const Matrix& des,const Matrix& tc, Matrix& cope,Mat
 
 float ols_dof(const Matrix& des){
   Matrix pdes = pinv(des);
-  Matrix R=Identity(des.Nrows())-des*pdes;
+  Matrix R=IdentityMatrix(des.Nrows())-des*pdes;
   return R.Trace();
 }