[Insight-users] Matrix multiplication bug

[Vincent Arsigny]

This is a multi-part message in MIME format.
--------------020201050609020906050907
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Content-Transfer-Encoding: 7bit

Hello everyone,

there is a bug in the following itk::Matrix operator :

vnl_matrix_fixed< T, NRows,
NColumns > operator * 
<http://www . itk . org/Insight/Doxygen/html/classitk_1_1Matrix . html#a4> 
(const vnl_matrix 
<http://www . itk . org/Insight/Doxygen/html/classvnl__matrix . html>< T > 
&matrix) const

You can find enclosed a very small test program + its CMakeList.txt that 
shows that the ITK * VNL matrix multiplication does work when the 
dimensions of the two matrices are identical but not in the general case 
like when multiplicating a 2*2 matrix by a 2*3 matrix.

Cheers,

Vincent Arsigny

--------------020201050609020906050907
Content-Type: text/plain;
 name="CMakeLists.txt"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="CMakeLists.txt"

#
# template for new module's CMakeList.txt
# this is an useless example only
#

FIND_PACKAGE(ITK)
IF (USE_ITK_FILE)
  INCLUDE(${USE_ITK_FILE})
ENDIF(USE_ITK_FILE)



#INCLUDE_DIRECTORIES(${MIPS_TRANSFORMATIONS_INCLUDE_DIR})



# Registration program with Levenberg-Marquardt alternating optimization,
# working in 2D only, specifically for 4 rigid components.



ADD_EXECUTABLE(TestSimple
	TestSimple.cxx
)

TARGET_LINK_LIBRARIES(TestSimple
   ITKIO
   ITKCommon
)
#ADD_EXECUTABLE(itkCommonTests
#	itkMatrixTest.cxx
#	itkCommonTests.cxx
#)

#TARGET_LINK_LIBRARIES(itkCommonTests
#   ITKIO
#   ITKCommon
#)


# Registration program with Levenberg-Marquardt alternating optimization,
# working in 2D only, specifically for 4 rigid components.
#
# This version enables fast registration, i.e. registration based on a
# limited subset of points, selected for their high gradient amplitude
# in the fixed image. This requires the storage in memory 
# of the gradient of the fixed image.


# fast registration is not ready yet. 27/11/03.


#(PolyRegistration2dLM-AltOpti-4C-FV
#	PolyRegistration2dLM-AltOpti-4C-FV.cxx
#)

#TARGET_LINK_LIBRARIES(PolyRegistration2dLM-AltOpti-4C-FV
#   ITKIO
#   ITKCommon
#   mipsPolyRigidRegistration
#)


# Parse the registration results for visualization purposes


#ADD_EXECUTABLE(
#  ReadPolyResults-4C
#  ReadPolyResults-4C.cxx
#  )

#TARGET_LINK_LIBRARIES(
#  ReadPolyResults-4C
#  ITKIO
#  ITKCommon
#  mipsPolyRigidRegistration
#)


# Parse the registration results for visualization purposes, special version
# to test fast registration


#ADD_EXECUTABLE(
#  ReadPolyResults-4C-FV
#  ReadPolyResults-4C-FV.cxx
#  )

#TARGET_LINK_LIBRARIES(
#  ReadPolyResults-4C-FV
#  ITKIO
#  ITKCommon
#  mipsPolyRigidRegistration
#)






--------------020201050609020906050907
Content-Type: text/plain;
 name="TestSimple.cxx"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="TestSimple.cxx"

#include "itkMatrix.h"
#include <iostream>
#include <fstream.h>



int main(int argc, char *argv[])
{
  std::cout << "\n\nSimple Test   starting now... !\n" << std::endl;
  
  try
    {
      

      itk::Matrix<double,2,2> m1;
      itk::Matrix<double,2,2> m2;
       
      for(int i=0;i<2;i++)
	for(int j=0;j<2;j++)
	  m1[i][j]=i+j;

      std::cout << "m1="  << std::endl;
      std::cout << m1 << std::endl;


      for(int i=0;i<2;i++)
	for(int j=0;j<2;j++)
	  m2[i][j]=i+j;

      std::cout << "m2=" << std::endl;
      std::cout << m2 << std::endl;


      std::cout << "VNL * VNL Multiplication result: " << std::endl;
      std::cout << m1.GetVnlMatrix()*m2.GetVnlMatrix() << std::endl;

      std::cout << "ITK * VNL Multiplication result: " << std::endl;
      std::cout << m1*m2.GetVnlMatrix() << std::endl;



      itk::Matrix<double,2,2> m3;
      itk::Matrix<double,2,3> m4;
       
      for(int i=0;i<2;i++)
	for(int j=0;j<2;j++)
	  m3[i][j]=i+j;

      std::cout << "m3="  << std::endl;
      std::cout << m3 << std::endl;


      for(int i=0;i<2;i++)
	for(int j=0;j<3;j++)
	  m4[i][j]=i+j;

      std::cout << "m4=" << std::endl;
      std::cout << m4 << std::endl;


      std::cout << "VNL * VNL Multiplication result: " << std::endl;
      std::cout << m3.GetVnlMatrix()*m4.GetVnlMatrix() << std::endl;

      std::cout << "ITK * VNL Multiplication result: " << std::endl;
      std::cout << m3*m4.GetVnlMatrix() << std::endl;
      


      
    } catch ( itk::ExceptionObject & e) {
	std::cout<<"Exception caught in test..."<<std::endl;
	std::cout<<e<<std::endl;
    }
      
  
  std::cout << "\nTest now finished.\n" << std::endl;

  return 0;
}

--------------020201050609020906050907--