[Insight-users] Matrix multiplication bug : now Fixed

[Luis Ibanez]

Hi Vincent,

Thanks for you detailed report.

You minimal code example made very easy for us to
locate an fix the bug.

The problem was in the retun type of the operator*().
It was using a vnl_matrix_fixed<> while it should be
using a vnl_matrix.

This has not been fixed.
http://www . itk . org/cgi-bin/cvsweb . cgi/Insight/Code/Common/itkMatrix . h?cvsroot=Insight&sortby=date

The problem was logged as Bug ID #412
http://www . itk . org/Bug/bug . php?op=show&bugid=412

Your nice test has now been added as part of itkMatrixTest.cxx
in Insight/Testing/Code/Common.

http://www . itk . org/cgi-bin/cvsweb . cgi/Insight/Testing/Code/Common/itkMatrixTest . cxx . diff?r1=1 . 16&r2=1 . 17&cvsroot=Insight&sortby=date


Note that you are also allowed to log bugs directly in the Bug
database. You can simply create an account at

            http://www . itk . org/Bug/

using your email as Id an making up your own password.
All users are welcome to log bug entries in this database.


Please let us know if you find any further problems.


Thanks,


    Luis


------------------------------------
Vincent Arsigny wrote:
> 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
> 
> 
> ------------------------------------------------------------------------
> 
> #
> # 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
> #)
> 
> 
> 
> 
> 
> 
> 
> ------------------------------------------------------------------------
> 
> #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;
> }