memory leak from mex file in R2012a?

Hi

At university I am using MATLAB version 7.14.0.739 (R2012a) on a Linux Ubuntu 32-bit system.

I am using code which calls a C++ mex file repeatedly. However, after a certain number of iterations I run out of memory.

As far as I can tell the C++ code, which was written by someone else, is being used correctly, and also I clear MATLAB variables as often as possible. But this still happens.

I don't know if I'll be able to upgrade my university version of MATLAB, so is there a memory leak known issue with mex files in R2012a?

Thanks

Matthew

Answer 1

[Philip Borghesani]

Most likely the problem is in the specific mex file. There are no systematic leaks with mex files in any recent version of MATLAB.

[Matthew Gillman]

Thanks.

The C++ code uses the mxCreateNumericArray() function for plhs a few times, and as soon as I have finished with it in the MATLAB code which calls the mex file I use the "clear" command. But the problem is still there...

[James Tursa]

Can you post your code, or a subset of your code, that produces the apparent problem? The only API function I know of that isn't garbage collected and can leak memory is mxArrayToString (it needs to be manually free'd with mxFree before the mex function returns to the caller).

SIDE NOTE to TMW: I think this is an oversight that should be corrected. Putting the address of the memory block returned by the mxArrayToString call on the garbage collection list should have been done in the first place IMO, and doing so won't break any code, so why not do it?

[Matthew Gillman]

Yes. I have run the memory profiler and it shows the MEX file (c_eoverlap) does not free any memory (in a MATLAB sense; there are delete[] statements in the C++). However, the MATLAB code which calls this MEX file, our_repeat.m, has a self memory which is negative and greater in magnitude by a greater amount than the MEX-allocated memory. That suggests the software is deleting the memory OK.

I do not know if the MEX file is a red herring when it comes to the memory leak. In a previous MATLAB file, I perform image transformations on a number of images and produce a dataset of transformed images for each original image. One transformation (blur) results in a dataset of 10 images for each original image; the other (light changes) yields datasets of size 14 images for each original. These images are then processed bythat MATLAB file. Why does this matter? Because, for each dataset, the MEX function is invoked for every member of that dataset; the memory allocation fails whilst processing the first image in the 22nd dataset. This is the case whether blur or light changes is used, i.e. it always fails on the same image despite the fact that a different number of calls have been made to the mex file, due to the differently-sized datasets.

But I post my code below for others to look at. I'm not including the driver code which runs the following.

(1.) find_repeatability_from_descriptors.m, written by me.

function [status] = find_repeatability_from_descriptors(category, ...\
    generated_files_location, generated_files_type, ...\
    path_to_results_directory, transformation, homography_directory)
% find_repeatability_from_descriptors.m
%
% example: category = 'abbey'
% generated_files_location = full path to the files produced by generate_files.m
% path_to_results_directory = directory where output file (stores repeatability values) will be produced.
% transformation - e.g. 'blur', 'light'
% homography_directory - directory where homography matrix is found
%
% A typical results file will be (e.g.) 'pgm_blur_abbey.csv'
%
% Original script: Matthew Gillman 01/06/2015.
% This version: 19/06/2015.
%
% This script is told the location of images. Each dataset consists of
% the original image and a number of images produced by
% applying a transformation of successive amounts on the original.
% These files will have been produced by generate_files.m
% That script also applied a number of point detectors to the images, 
% and then ran a point descriptor binary on these. 
% It is assumed that the descriptor files are in the same directory as the image files.
% This script takes the image file datasets and descriptor files and finds the repeatability
% between each original image and one of its transformed images.  
%
% "A software that takes images of a scene category as input, then apply
% the desired transformation on each image of the scene category to make a
% number of datasets, and then run all feature detectors on these datasets
% to generate repeatability values, and finally repeatability curves
% (as in the thesis)."
%

%
% TODO
% Add details here of params etc.
% Add brief pause() to loop to better enable CTRL-C handling.

      generated_files_location    % used for jpg files, detector files, etc.
      generated_files_type
      path_to_results_directory
      transformation % 'blur', 'compress', 'light' 
      status = 0;

      switch generated_files_type
          case {'jpg', 'JPG', 'jpeg', 'JPEG'}
              generated_files_ext = '.jpg';
          case {'pgm', 'PGM'}
              generated_files_ext = '.pgm';
          otherwise
              fprintf('\nError: I cannot handle file type: %s\n', generated_files_ext);
              return;
      end;

      % Store step values required for our transformation. The precise
      % values will depend on the transformation chosen.
      % Seeing as the following are used in both generate_files_1.m and this
      % script, maybe it would be better to pass in the values of the 
      % transformations required to both scripts, rather than duplicating
      % code here.
      transvals = [];    

      switch transformation
              case 'blur'
                  transvals = 0.5:0.5:4.5;
              case 'light'
                  transvals = [0.9 0.8 0.7 0.6 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1];
              otherwise
                  fprintf('\nSorry, I cannot handle %s\n', transformation);
                  return;
      end

      %reference_image_files = dir([generated_files_location, '/sun*.jpg']);
      reference_image_files = dir([generated_files_location, '/sun*', generated_files_ext]);
      %reference_image_files
      %return;
      len = length(reference_image_files); % No. of datasets.

      % Store the files produced by the transformations. Untransformed image
      % goes in index 1 and transformed ones subsequently.
      % Have to use cell arrays 
      % http://stackoverflow.com/questions/2288899/appending-string-to-matlab-array
      % Each row (dataset) of all_image_files stores 1 original image file and its
      % transformed counterparts.
      % Pre-allocate memory for speed.
      all_image_files = cell(len, length(transvals) + 1);
      %all_image_files
      % index = 1;
      %y = reference_image_files(1).name
      %z = reference_image_files(2).name

      % Have to use a loop because I can't get the quick way to work.
      for i=1:len
         all_image_files(i,1) = {reference_image_files(i).name};
      end

      %all_image_files(1,:)  = {reference_image_files.name};
  %     all_image_files(1,1)
  %     all_image_files(2,1)    
  %     all_image_files(3,1)
  %     reference_image_files
      %all_image_files
      % y = all_image_files{3,1};
      %y.name
      %all_image_files(1,1).name
     % all_image_files.name(1,1)

      clear reference_image_files;

      % Now add the transformed jpg files to the cell array.
      for i = 1:len
          index = 1;
          reference = all_image_files{i,1};
          %reference
          %reference.name

          for t = transvals(1:end)
              index = index + 1;
              switch transformation
                  case 'blur'
                      value = t; % e.g. 0.5, 1, 1.5...
                  case 'light'
                      value = sprintf('%3.2f', 1-t); % e.g. 0.10
                  otherwise
                      printf('\nERROR: I cannot handle %s!\n', transformation);
              end

              %filename = sprintf('%s_%s', num2str(value), char({reference})); %[num2str(value), '_', reference]
              %filename = [reference];
              %reference(i,1).name
              filename = [num2str(value) '_' reference];
              filename
              all_image_files(i, index) = {filename}; % e.g. '4_sun_aaalbzqrimafwbiv.jpg'

          end % t

      end % i

      % Create output .csv file for writing. 
      %
      % Thus a typical file name will be (say) pgm_blur_abbey.csv or perhaps
      % jpg_blur_abbey.csv. This stops files with the same name (e.g. just 
      % 'blur_abbey.csv') overwriting each other (i.e. if you have both jpg
      % and pgm versions, required for the different detector types).
      % 
      output_file_name = [generated_files_type, '_', transformation, '_', category, '.csv'];
      outfile = fullfile(path_to_results_directory, output_file_name);
      fileID=fopen(outfile, 'wt');
      if ( fileID < 3 )
          fprintf('\nERROR: could not open %s for writing; fileID = %d\n', outfile, fileID);
          return;
      end

      write_output_file_header(fileID, transformation, transvals);

      % Find the repeatability for the different detectors.
      % Note that it is OK to pass the large cell array all_image_files to the
      % function. From
      % http://uk.mathworks.com/matlabcentral/ans return;wers/152-can-matlab-pass-by-reference:
      % "If instead you are attempting to use pass-by-reference to avoid
      % unnecessary copying of data into the workspace of the function you're
      % calling, you should be aware that MATLAB uses a system commonly called
      % 'copy-on-write' to avoid making a copy of the input argument inside
      % the function workspace until or unless you modify the input argument.
      % If you do not modify the input argument, MATLAB will avoid making a
      % copy."

      % Not all image file types are compatible with all detector binaries,
      % so we only slurp up the files we need generated by these.

      switch generated_files_ext
         case '.jpg'
              get_repeatability('MSER', fileID, generated_files_location, generated_files_ext, all_image_files, homography_directory);
              get_repeatability('EBR', fileID, generated_files_location, generated_files_ext, all_image_files, homography_directory);
              get_repeatability('IBR', fileID, generated_files_location, generated_files_ext, all_image_files, homography_directory);
         case '.pgm'
              get_repeatability('Harris-Affine', fileID, generated_files_location, generated_files_ext, all_image_files, homography_directory);
              get_repeatability('Hessian-Affine', fileID, generated_files_location, generated_files_ext, all_image_files, homography_directory);
              %get_repeatability('Harris', fileID, generated_files_location, generated_files_ext, all_image_files, homography_directory);
              %get_repeatability('Hessian-Laplace', fileID, generated_files_location, generated_files_ext, all_image_files, homography_directory);
              %get_repeatability('Harris-Laplace', fileID, generated_files_location, generated_files_ext, all_image_files, homography_directory);
              % surf
         otherwise
              fprintf('\nI should never reach here! File ext is: %s\n', generated_files_ext);
              return;
     end

     fclose(fileID);

     %TODO would be better for both scripts to take in the transformation values as a parameter.

end % find_repeatability()

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% write_output_file_header()
%
% Given a transformation transf (e.g. 'blur') and set of 
% values (tvals) for this transformation, write the header
% of the output file.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function write_output_file_header(fileID, transf, tvals)

      fprintf(fileID, '%s', '');
      for value = tvals(1:end)

          string = [',Repeatability for ', num2str(value), ' ', transf]; 

          %value
          %string

          fprintf(fileID, '%s', string);

      end

      fprintf(fileID, '\n');
      %flushoutput(fileID);%function too new

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% get_repeatability()
%
% For the given detector, load in the descriptor files and 
% compute the repeatability. Write
% the results to the output file.
%
% detector =
% 'Harris'|'Harris-Laplace'|'Hessian-Laplace'|'Harris-Affine'|
% 'Hessian-Affine'|'IBR'|'EBR'|'MSER'
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [status] = get_repeatability(detector, fileID, generated_files_location, generated_files_ext, all_image_files, homog_dir)

      fprintf(fileID, '%s', detector);
      fprintf(fileID, '\n\n');

      sz = size(all_image_files);
      nrows = sz(1)
      num_per_dataset = sz(2)

      % TODO: the next line is only valid in cases like blurring, lighting
      % and jpeg compression, i.e. when there is no camera movement.
      H = fullfile(homog_dir, 'Identity'); % Homography matrix.

      % Loop over every file (dataset = original file and its transformed
      % ones).
      % Produce 1 line per dataset in output file.
      % NB some files will be missing if the detectors failed to find any
      % points in the corresponding transformed image file. We give those
      % a repeatability of 0.

      for i=1:nrows

          dataset = all_image_files(i,:); % Reference file + its transformed files.
          empty = cellfun('isempty', dataset);

          % Remove '.jpg' or '.pgm' from filenames. e.g. 'myfile.jpg'->'myfile'
          stubs = cell(1, num_per_dataset); 
          %stubs = strrep(dataset, '.jpg', ''); 
          % get error with above if empty cells at end of array.
          for j=1:num_per_dataset
             if ~empty(j)
                 mycell = dataset(j);
                 stubs{j} = strrep(mycell, generated_files_ext, ''); 
             end
          end
          %stubs{1}
          %stubs{10}

          % Convert stubs cell array to matrix and pull off first item
          % (i.e. the reference file without its .jpg or .pgm extension)
          % From http://uk.mathworks.com/matlabcentral/answers/20281-cell-to-matrix
          %stubmatrix = cat(1, stubs);
          %reference = stubmatrix(1,:);
          %stubmatrix = cell2mat(stubs)
          %reference = stubmatrix(1)
          reference = cell2mat(stubs{1});
          %reference = stubs{1};
          %reference
          %return;
          fprintf(fileID, '%s', reference);
          points_ext = []; % Points file extension.

          switch detector % TODO: Add salient etc.
              case 'Harris' 
                  points_ext = '.har';
              case 'Harris-Laplace'
                  points_ext = '.harlap';
              case 'Hessian-Laplace'
                  points_ext = '.heslap';
              case 'Harris-Affine'
                  points_ext = '.haraff';
              case 'Hessian-Affine'
                  points_ext = '.hesaff';
              case 'EBR'
                  points_ext = '.ebr';
              case 'IBR'
                  points_ext = '.ibr';
              case 'MSER'
                  points_ext = '.mser';
              otherwise
                  fprintf('\nERROR: get_repeatability() cannot handle detector %s\n', detector);
                  return;
          end % switch

          points = strcat(stubs, points_ext);  % e.g. 'myfile.har'
          %sifts = strcat(stubs, 'haraff.sift');
          %sifts = strcat(points, '.sift');     % e.g. 'myfile.har.sift'

          full_images{1, num_per_dataset} = []; % Full paths to .jpg files
          full_points{1, num_per_dataset} = []; % Full path to (e.g.) .har files

          for j = 1:num_per_dataset
              %stub = stubs(j);   %stubs(1)
              %stubs(3)

              %fullstub = fullfile(temp_files_location, stub);

              %img = cell2mat(dataset(j))
              %imgpath = fullfile(temp_files_location, img)
              %full_images{j} = imgpath
              if ~empty(j)
                  full_images{j} = fullfile(generated_files_location, cell2mat(dataset(j)));
                  full_points{j} = fullfile(generated_files_location, cell2mat(points{j}));
              end % if ~empty(j)

          end % for j = 1:num_per_dataset

          % We have now found the descriptors, so obtain repeatability and
          % other info.matches
          % Write the number of features found in the common region of image 1
          % and image 2, and the number of repeated features in image 2, to the
          % file.
          % [erro,repeat,corresp, match_score,matches, twi,s1,s2,common1,common2]=our_repeat(file1,file2,H,imf1,imf2,1);
          % Note that even though images transformed by increasing amounts will
          % all be compared with the original image, "common1" and "common2" will
          % usually be different every time the comparisons are made. 

          %full_images(1)
          %full_images{4}

          for j=2:num_per_dataset

              if ~empty(j)

                  if exist(full_points{j}, 'file') == 0
                      repeat_results = 0;
                  else

                      [erro,repeat,corresp, match_score,matches, twi,s1,s2,common1,common2] ...\
                          =our_repeat(full_points{1}, full_points{j}, ...\
                              H, full_images{1}, full_images{j}, 1);    

                      %corresp_results = corresp(4);
                      repeat_results = repeat(4);
                      clear erro repeat corresp match_score matches twi s1 s2 common1 common2;

                  end

                  % TODO: waiting for Shoaib to confirm these are what's required.
                  %fprintf(fileID, ',%d,%d,%d,%d',repeat_results, common1, common2, matches);
                  %fprintf(fileID, ',%d', repeat_results); e.g. 7.447584e+01
                  fprintf(fileID, ',%.2f', repeat_results); % e.g. 74.48                          

              end % if

              clear repeat_results;

          end % for j

          fprintf(fileID, '\n');
          %flushoutput(fileID);

          clear dataset full_images full_points;

      end % for i=1:nrows

      fprintf(fileID, '\n');

end % get_repeatability()
****************************************************************

    (2.) our_repeat.m, written by someone else:

function [erro,repeat,corresp, match_score,matches, twi,s1,s2, common1, common2]=our_repeat(file1,file2,Hom,imf1,imf2, common_part)
%
%
%Computes repeatability and overlap score between two lists of features
%detected in two images.
%   [erro,repeat,corresp,matched,matchedp]=repeatability('file1','file2','H4to1','imf1','imf2','-ko');
%
%IN: 
%    file1 - file 1 with detected features         
%    file2 - file 2 with detected features 
%    H - file with 3x3 Homography matrix from image 1 to 2, x2=H*x1
%        Assumes that image coordiantes are 0..width.
%    imf1 - image file  1
%    imf2 - image file  2
%    common_part - flag should be set to 1 for repeatability and 0 for descriptor performance
%
%OUT :    erro - overlap %
%         repeat - repeatability %
%         corresp - nb of correspondences
%         match_score  - matching score
%         matches - number of correct nearest neighbour matches
%         twi - matrix with overlap errors<50\%
%
%  region file format :
%--------------------
%descriptor_size  
%nbr_of_regions
%x1 y1 a1 b1 c1 d1 d2 d3 ...
%x2 y2 a2 b2 c2 d1 d2 d3 ...
%....
%....
%---------------------
%x, y - center coordinates
%a, b, c - ellipse parameters ax^2+2bxy+cy^2=1
%d1 d2 d3 ... - descriptor invariants
%if descriptor_size<=1 the descriptor is ignored

fprintf(1,'Reading and sorting the regions...\n');

[f1 s1 dimdesc1]=loadFeatures(file1);
[f2 s2 dimdesc2]=loadFeatures(file2);

H=load(Hom);

fprintf(1,'nb of regions in file1 %d - descriptor dimension %d.\n',s1,dimdesc1);
fprintf(1,'nb of regions in file2 %d - descriptor dimension %d.\n',s2,dimdesc2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% my_no_features_img1 = s1;
% my_no_features_img2 = s2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Added by Matthew Gillman 03/06/2015:
% common1 = no. of regions in common part of image1
% common2 = ditto image2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if size(f1,1)==5 & size(f1,1)==size(f2,1) 
fprintf(1,'%s looks like file with affine regions...\n',file1);
  if  size(f1,1)~= 5 | size(f1,1) ~= 5
    error('Wrong ascii format of %s or %s files.',file1,file2);
  end
elseif dimdesc1>1 & dimdesc1==dimdesc2 
  fprintf(1,'%s, %s look like files with descriptors...\n',file1,file2);
elseif s1 == 0 | s2 == 0 % Added by Matthew Gillman
  fprintf('\nZero descriptors found in %s and/or %s; repeatability will be set to 0\n', file1, file2);
  erro = [10:10:60];
  corresp = zeros(1,6);
  repeat = zeros(1,6);
  matches = 0;
  match_score = 0;
  twi = [];
  common1 = 0;
  common2 = 0;
  return;
else
    error('Different descriptor dimension in %s or %s files.',file1,file2);
end

dimdesc=dimdesc1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% allocate vectors for the features
feat1=zeros(9+dimdesc,s1);
feat2=zeros(9+dimdesc,s2);
feat1(1:5,1:s1)=f1(1:5,1:s1);
feat2(1:5,1:s2)=f2(1:5,1:s2);
if size(f1,1)>1
feat1(10:9+dimdesc,1:s1)=f1(6:5+dimdesc,1:s1);
feat2(10:9+dimdesc,1:s2)=f2(6:5+dimdesc,1:s2);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%project regions from one image to the other
HI=H(:, 1:3);
H=inv(HI);
fprintf(1,'Projecting 1 to 2...');
[feat1 feat1t scales1]=project_regions(feat1',HI);
fprintf(1,'and 2 to 1...\n');
[feat2 feat2t scales2]=project_regions(feat2',H);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if common_part==1
fprintf(1,'Removing features from outside of the common image part...\n');
im1=imread(imf1);
im2=imread(imf2);
im1x=size(im1);
im1y=im1x(1);
im1x=im1x(2);
im2x=size(im2);
im2y=im2x(1);
im2x=im2x(2);
ind=find((feat1(:,1)+feat1(:,8))<im1x & (feat1(:,1)-feat1(:,8))>0 & (feat1(:,2)+feat1(:,9))<im1y & (feat1(:,2)-feat1(:,9))>0);
feat1=feat1(ind,:);
feat1t=feat1t(ind,:);
ind=find((feat1t(:,1)+feat1t(:,8))<im2x & (feat1t(:,1)-feat1t(:,8))>0 & (feat1t(:,2)+feat1t(:,9))<im2y & (feat1t(:,2)-feat1t(:,9))>0);
feat1=feat1(ind,:);
feat1t=feat1t(ind,:);
scales1=scales1(ind);

ind=find((feat2(:,1)+feat2(:,8))<im2x & (feat2(:,1)-feat2(:,8))>0 & (feat2(:,2)+feat2(:,9))<im2y & (feat2(:,2)-feat2(:,9))>0);
feat2t=feat2t(ind,:);
feat2=feat2(ind,:);
ind=find((feat2t(:,1)+feat2t(:,8))<im1x & (feat2t(:,1)-feat2t(:,8))>0 & (feat2t(:,2)+feat2t(:,9))<im1y & (feat2t(:,2)-feat2t(:,9))>0);
feat2t=feat2t(ind,:);
feat2=feat2(ind,:);
scales2=scales2(ind);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

common1 = size(feat1,1);
common2 = size(feat2t,1);
fprintf(1,'nb of regions in common part in image1 %d.\n', common1);
fprintf(1,'nb of regions in common part in image2 %d.\n', common2);
%fprintf(1,'nb of regions in common part in image1 %d.\n',size(feat1,1));
%fprintf(1,'nb of regions in common part in image2 %d.\n',size(feat2t,1));
end 

%sf=min([size(feat1,1) size(feat2t,1)]);
sf= size(feat1,1);

feat1=feat1';
feat1t=feat1t';
feat2t=feat2t';
feat2=feat2';

fprintf(1,'Computing overlap error & selecting one-to-one correspondences: ');
tic;
%c_eoverlap is a C implementation to compute the overlap error.
[wout, twout, dout, tdout]=c_eoverlap(feat1,feat2t,common_part);
t=toc;
fprintf(1,' %.1f sec.\n',t);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

erro=[10:10:60];
corresp=zeros(1,6);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%compute the number of correspondences
for i=1:6,
wi=(wout<erro(i));
corresp(i)=sum(sum(wi));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
repeat=100*corresp/sf;

fprintf(1,'\noverlap error: ');
fprintf(1,'%.1f ',erro);
fprintf(1,'\nrepeatability: ');
fprintf(1,'%.1f ',repeat);
fprintf(1,'\nnb of correspondences: ');
fprintf(1,'%.0f ',corresp);
fprintf(1,'\n');

match_overlap=40;
if common_part==0
match_overlap=50;
end
fprintf(1,'Matching with the descriptor for the overlap error < %d%\n',match_overlap);
if common_part==0
twi=(twout<match_overlap);
else
twi=(wout<match_overlap);
end

dx=(dout<10000).*(twi);
matches=sum(sum(dx));
match_score=100*matches/sf;

fprintf(1,'\nMatching score  %0.1f, nb of correct matches %.1f.\n',match_score,matches);

% Added by Matthew Gillman
clear H HI feat1t feat2t scales1 scales2 feat1 feat2 f1 f2 dimdesc1 dimdesc2 ...\
 wout twout dout tdout dx match_overlap wi im1x im1y im2x im2y ind;

%end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [feat,featp,scales]=project_regions(feat,H)

s=size(feat);
s1=s(1);

featp=feat;
scales=zeros(1,s1);

for c1=1:s1,%%%%%%%%%%%%%%%%%
%feat(c1,3:5)=(1/25)*feat(c1,3:5);
Mi1=[feat(c1,3) feat(c1,4);feat(c1,4) feat(c1,5)];

%compute affine transformation
[v1 e1]=eig(Mi1);
d1=(1/sqrt(e1(1))); 
d2=(1/sqrt(e1(4))); 
sc1=sqrt(d1*d2);
feat(c1,6)=d1;
feat(c1,7)=d2; 
scales(c1)=sqrt(feat(c1,6)*feat(c1,7));

%bounding box
feat(c1,8) = sqrt(feat(c1,5)/(feat(c1,3)*feat(c1,5) - feat(c1,4)^2));
feat(c1,9) = sqrt(feat(c1,3)/(feat(c1,3)*feat(c1,5) - feat(c1,4)^2));

Aff=getAff(feat(c1,1),feat(c1,2),sc1, H);

%project to image 2
l1=[feat(c1,1),feat(c1,2),1];function [erro,repeat,corresp, match_score,matches, twi,s1,s2, common1, common2]=our_repeat(file1,file2,Hom,imf1,imf2, common_part)
%
%
%Computes repeatability and overlap score between two lists of features
%detected in two images.
%   [erro,repeat,corresp,matched,matchedp]=repeatability('file1','file2','H4to1','imf1','imf2','-ko');
%
%IN: 
%    file1 - file 1 with detected features         
%    file2 - file 2 with detected features 
%    H - file with 3x3 Homography matrix from image 1 to 2, x2=H*x1
%        Assumes that image coordiantes are 0..width.
%    imf1 - image file  1
%    imf2 - image file  2
%    common_part - flag should be set to 1 for repeatability and 0 for descriptor performance
%
%OUT :    erro - overlap %
%         repeat - repeatability %
%         corresp - nb of correspondences
%         match_score  - matching score
%         matches - number of correct nearest neighbour matches
%         twi - matrix with overlap errors<50\%
%
%  region file format :
%--------------------
%descriptor_size  
%nbr_of_regions
%x1 y1 a1 b1 c1 d1 d2 d3 ...
%x2 y2 a2 b2 c2 d1 d2 d3 ...
%....
%....
%---------------------
%x, y - center coordinates
%a, b, c - ellipse parameters ax^2+2bxy+cy^2=1
%d1 d2 d3 ... - descriptor invariants
%if descriptor_size<=1 the descriptor is ignored

fprintf(1,'Reading and sorting the regions...\n');

[f1 s1 dimdesc1]=loadFeatures(file1);
[f2 s2 dimdesc2]=loadFeatures(file2);

H=load(Hom);

fprintf(1,'nb of regions in file1 %d - descriptor dimension %d.\n',s1,dimdesc1);
fprintf(1,'nb of regions in file2 %d - descriptor dimension %d.\n',s2,dimdesc2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% my_no_features_img1 = s1;
% my_no_features_img2 = s2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Added by Matthew Gillman 03/06/2015:
% common1 = no. of regions in common part of image1
% common2 = ditto image2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if size(f1,1)==5 & size(f1,1)==size(f2,1) 
fprintf(1,'%s looks like file with affine regions...\n',file1);
  if  size(f1,1)~= 5 | size(f1,1) ~= 5
    error('Wrong ascii format of %s or %s files.',file1,file2);
  end
elseif dimdesc1>1 & dimdesc1==dimdesc2 
  fprintf(1,'%s, %s look like files with descriptors...\n',file1,file2);
elseif s1 == 0 | s2 == 0 % Added by Matthew Gillman
  fprintf('\nZero descriptors found in %s and/or %s; repeatability will be set to 0\n', file1, file2);
  erro = [10:10:60];
  corresp = zeros(1,6);
  repeat = zeros(1,6);
  matches = 0;
  match_score = 0;
  twi = [];
  common1 = 0;
  common2 = 0;
  return;
else
    error('Different descriptor dimension in %s or %s files.',file1,file2);
end

dimdesc=dimdesc1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% allocate vectors for the features
feat1=zeros(9+dimdesc,s1);
feat2=zeros(9+dimdesc,s2);
feat1(1:5,1:s1)=f1(1:5,1:s1);
feat2(1:5,1:s2)=f2(1:5,1:s2);
if size(f1,1)>1
feat1(10:9+dimdesc,1:s1)=f1(6:5+dimdesc,1:s1);
feat2(10:9+dimdesc,1:s2)=f2(6:5+dimdesc,1:s2);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%project regions from one image to the other
HI=H(:, 1:3);
H=inv(HI);
fprintf(1,'Projecting 1 to 2...');
[feat1 feat1t scales1]=project_regions(feat1',HI);
fprintf(1,'and 2 to 1...\n');
[feat2 feat2t scales2]=project_regions(feat2',H);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if common_part==1
fprintf(1,'Removing features from outside of the common image part...\n');
im1=imread(imf1);
im2=imread(imf2);
im1x=size(im1);
im1y=im1x(1);
im1x=im1x(2);
im2x=size(im2);
im2y=im2x(1);
im2x=im2x(2);
ind=find((feat1(:,1)+feat1(:,8))<im1x & (feat1(:,1)-feat1(:,8))>0 & (feat1(:,2)+feat1(:,9))<im1y & (feat1(:,2)-feat1(:,9))>0);
feat1=feat1(ind,:);
feat1t=feat1t(ind,:);
ind=find((feat1t(:,1)+feat1t(:,8))<im2x & (feat1t(:,1)-feat1t(:,8))>0 & (feat1t(:,2)+feat1t(:,9))<im2y & (feat1t(:,2)-feat1t(:,9))>0);
feat1=feat1(ind,:);
feat1t=feat1t(ind,:);
scales1=scales1(ind);

ind=find((feat2(:,1)+feat2(:,8))<im2x & (feat2(:,1)-feat2(:,8))>0 & (feat2(:,2)+feat2(:,9))<im2y & (feat2(:,2)-feat2(:,9))>0);
feat2t=feat2t(ind,:);
feat2=feat2(ind,:);
ind=find((feat2t(:,1)+feat2t(:,8))<im1x & (feat2t(:,1)-feat2t(:,8))>0 & (feat2t(:,2)+feat2t(:,9))<im1y & (feat2t(:,2)-feat2t(:,9))>0);
feat2t=feat2t(ind,:);
feat2=feat2(ind,:);
scales2=scales2(ind);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

common1 = size(feat1,1);
common2 = size(feat2t,1);
fprintf(1,'nb of regions in common part in image1 %d.\n', common1);
fprintf(1,'nb of regions in common part in image2 %d.\n', common2);
%fprintf(1,'nb of regions in common part in image1 %d.\n',size(feat1,1));
%fprintf(1,'nb of regions in common part in image2 %d.\n',size(feat2t,1));
end 

%sf=min([size(feat1,1) size(feat2t,1)]);
sf= size(feat1,1);

feat1=feat1';
feat1t=feat1t';
feat2t=feat2t';
feat2=feat2';

fprintf(1,'Computing overlap error & selecting one-to-one correspondences: ');
tic;
%c_eoverlap is a C implementation to compute the overlap error.
[wout, twout, dout, tdout]=c_eoverlap(feat1,feat2t,common_part);
t=toc;
fprintf(1,' %.1f sec.\n',t);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

erro=[10:10:60];
corresp=zeros(1,6);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%compute the number of correspondences
for i=1:6,
wi=(wout<erro(i));
corresp(i)=sum(sum(wi));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
repeat=100*corresp/sf;

fprintf(1,'\noverlap error: ');
fprintf(1,'%.1f ',erro);
fprintf(1,'\nrepeatability: ');
fprintf(1,'%.1f ',repeat);
fprintf(1,'\nnb of correspondences: ');
fprintf(1,'%.0f ',corresp);
fprintf(1,'\n');

match_overlap=40;
if common_part==0
match_overlap=50;
end
fprintf(1,'Matching with the descriptor for the overlap error < %d%\n',match_overlap);
if common_part==0
twi=(twout<match_overlap);
else
twi=(wout<match_overlap);
end

dx=(dout<10000).*(twi);
matches=sum(sum(dx));
match_score=100*matches/sf;

fprintf(1,'\nMatching score  %0.1f, nb of correct matches %.1f.\n',match_score,matches);

% Added by Matthew Gillman
clear H HI feat1t feat2t scales1 scales2 feat1 feat2 f1 f2 dimdesc1 dimdesc2 ...\
 wout twout dout tdout dx match_overlap wi im1x im1y im2x im2y ind;

%end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [feat,featp,scales]=project_regions(feat,H)

s=size(feat);
s1=s(1);

featp=feat;
scales=zeros(1,s1);

for c1=1:s1,%%%%%%%%%%%%%%%%%
%feat(c1,3:5)=(1/25)*feat(c1,3:5);
Mi1=[feat(c1,3) feat(c1,4);feat(c1,4) feat(c1,5)];

%compute affine transformation
[v1 e1]=eig(Mi1);
d1=(1/sqrt(e1(1))); 
d2=(1/sqrt(e1(4))); 
sc1=sqrt(d1*d2);
feat(c1,6)=d1;
feat(c1,7)=d2; 
scales(c1)=sqrt(feat(c1,6)*feat(c1,7));

%bounding box
feat(c1,8) = sqrt(feat(c1,5)/(feat(c1,3)*feat(c1,5) - feat(c1,4)^2));
feat(c1,9) = sqrt(feat(c1,3)/(feat(c1,3)*feat(c1,5) - feat(c1,4)^2));

Aff=getAff(feat(c1,1),feat(c1,2),sc1, H);

%project to image 2
l1=[feat(c1,1),feat(c1,2),1];
l1_2=H*l1';
l1_2=l1_2/l1_2(3);
featp(c1,1)=l1_2(1);
featp(c1,2)=l1_2(2);
BMB=inv(Aff*inv(Mi1)*Aff');
[v1 e1]=eig(BMB);
featp(c1,6)=(1/sqrt(e1(1)));
featp(c1,7)=(1/sqrt(e1(4))); 
featp(c1,3:5)=[BMB(1) BMB(2) BMB(4)];
%bounding box in image 2
featp(c1,8) = sqrt(featp(c1,5)/(featp(c1,3)*featp(c1,5) - featp(c1,4)^2));
featp(c1,9) = sqrt(featp(c1,3)/(featp(c1,3)*featp(c1,5) - featp(c1,4)^2));
% Added by Matthew
clear Aff BMB Mi1;
end
%end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Aff=getAff(x,y,sc,H)
h11=H(1);
h12=H(4);
h13=H(7);
h21=H(2);
h22=H(5);
h23=H(8);
h31=H(3);
h32=H(6);
h33=H(9);
fxdx=h11/(h31*x + h32*y +h33) - (h11*x + h12*y +h13)*h31/(h31*x + h32*y +h33)^2;
fxdy=h12/(h31*x + h32*y +h33) - (h11*x + h12*y +h13)*h32/(h31*x + h32*y +h33)^2;

fydx=h21/(h31*x + h32*y +h33) - (h21*x + h22*y +h23)*h31/(h31*x + h32*y +h33)^2;
fydy=h22/(h31*x + h32*y +h33) - (h21*x + h22*y +h23)*h32/(h31*x + h32*y +h33)^2;

            Aff=[fxdx fxdy;fydx fydy];
  %end

function [feat nb dim]=loadFeatures(file)
fid = fopen(file, 'r');
dim=fscanf(fid, '%f',1);
if dim==1
dim=0;
end
nb=fscanf(fid, '%d',1);
feat = fscanf(fid, '%f', [5+dim, inf]);
fclose(fid);
%end

l1_2=H*l1';
l1_2=l1_2/l1_2(3);
featp(c1,1)=l1_2(1);
featp(c1,2)=l1_2(2);
BMB=inv(Aff*inv(Mi1)*Aff');
[v1 e1]=eig(BMB);
featp(c1,6)=(1/sqrt(e1(1)));
featp(c1,7)=(1/sqrt(e1(4))); 
featp(c1,3:5)=[BMB(1) BMB(2) BMB(4)];
%bounding box in image 2
featp(c1,8) = sqrt(featp(c1,5)/(featp(c1,3)*featp(c1,5) - featp(c1,4)^2));
featp(c1,9) = sqrt(featp(c1,3)/(featp(c1,3)*featp(c1,5) - featp(c1,4)^2));
% Added by Matthew
clear Aff BMB Mi1;
end
%end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Aff=getAff(x,y,sc,H)
h11=H(1);
h12=H(4);
h13=H(7);
h21=H(2);
h22=H(5);
h23=H(8);
h31=H(3);
h32=H(6);
h33=H(9);
fxdx=h11/(h31*x + h32*y +h33) - (h11*x + h12*y +h13)*h31/(h31*x + h32*y +h33)^2;
fxdy=h12/(h31*x + h32*y +h33) - (h11*x + h12*y +h13)*h32/(h31*x + h32*y +h33)^2;

fydx=h21/(h31*x + h32*y +h33) - (h21*x + h22*y +h23)*h31/(h31*x + h32*y +h33)^2;
fydy=h22/(h31*x + h32*y +h33) - (h21*x + h22*y +h23)*h32/(h31*x + h32*y +h33)^2;

            Aff=[fxdx fxdy;fydx fydy];
  %end

function [feat nb dim]=loadFeatures(file)
fid = fopen(file, 'r');
dim=fscanf(fid, '%f',1);
if dim==1
dim=0;
end
nb=fscanf(fid, '%d',1);
feat = fscanf(fid, '%f', [5+dim, inf]);
fclose(fid);
%end
***************************************************************

(3.) C++ MEX code, written by someone else:

#include <mex.h>
#include <math.h>

// This mex function implements the equivalent of:
//
//   function y = example_mat(x)
//   y = 2 * x + 1;
//   return

// [y, z] = foo(a, b, c)
//
//  nrhs = 3
//  a : prhs[0]
//  b : prhs[1]
//  c : prhs[2]
//
//  nlhs = 2
//  y : plhs[0]
//  z : plhs[1]

void mexFunction(int nlhs,       mxArray *plhs[],
                 int nrhs, const mxArray *prhs[])
{
  // must have single double array input.
  if (nrhs != 3) { printf("in nrhs != 3\n"); return; }
  if (mxGetClassID(prhs[0]) != mxDOUBLE_CLASS) { printf("input must be a double array\n"); return; }
  if (mxGetClassID(prhs[1]) != mxDOUBLE_CLASS) { printf("input must be a double array\n"); return; }

    // must have single output.
    if (nlhs != 4) { printf("out nlhs != 4\n"); return; }

    // get size of x.
    int const num_dims1 = mxGetNumberOfDimensions(prhs[0]);
    int const *dims1 = mxGetDimensions(prhs[0]);
    int const num_dims2 = mxGetNumberOfDimensions(prhs[1]);
    int const *dims2 = mxGetDimensions(prhs[1]);
    int const num_dims3 = mxGetNumberOfDimensions(prhs[2]);
    int const *dims3 = mxGetDimensions(prhs[2]);
    if(dims1[0]<9 || dims2[0]<9 || dims1[0]!=dims2[0] || dims3[0]!=1){ printf("dims != 9\n"); return; }

    // create new array of the same size as x.

    mwSize /*int*/ *odim = new mwSize /*int*/ [2];
    odim[0]=dims2[1];
    odim[1]=dims1[1];
    plhs[0] = mxCreateNumericArray(2, odim, mxDOUBLE_CLASS, mxREAL);
    plhs[1] = mxCreateNumericArray(2, odim, mxDOUBLE_CLASS, mxREAL);
    plhs[2] = mxCreateNumericArray(2, odim, mxDOUBLE_CLASS, mxREAL);
    plhs[3] = mxCreateNumericArray(2, odim, mxDOUBLE_CLASS, mxREAL);

    // get pointers to beginning of x and y.
    double const *feat1 = (double *)mxGetData(prhs[0]);
    double const *feat2 = (double *)mxGetData(prhs[1]);
    double const *flag = (double *)mxGetData(prhs[2]);
  // Following changes made by Matthew Gillman after seeing post on Mathworks
  // website in "Creating Potential Memory Leaks" at
  // uk.mathworks.com/help/matlab/matlab_external/memory-management-issues.html
  // unfortunately this makes no difference
   // double       *over_out = (double *)mxGetData(plhs[0]);
    double *over_out = mxGetPr(plhs[0]);  
    //double       *mover_out = (double *)mxGetData(plhs[1]);
  double *mover_out = mxGetPr(plhs[1]);
    //double       *desc_out = (double *)mxGetData(plhs[2]);
  double *desc_out = mxGetPr(plhs[2]);
    //double       *mdesc_out = (double *)mxGetData(plhs[3]);
  double *mdesc_out = mxGetPr(plhs[3]);

    float *feat1a = new float[9];
    float *feat2a = new float[9];
    float *tdesc_out = new float[dims2[1]*dims1[1]];
    float *tover_out = new float[dims2[1]*dims1[1]];

    int common_part=(int)flag[0];

     for(int j=0;j<dims1[1];j++){    
      for (int i=0; i<dims2[1]; i++){
        over_out[j*dims2[1]+i]=100.0;
        desc_out[j*dims2[1]+i]=1000000.0;
      }
     } 

     // printf("%f %f\n",flag[0],flag[1]);
    // total number of elements in arrays.
    /*int total = 1;
    for (int i=0; i<num_dims1; ++i){
      total *= dims1[i];  
      printf("feat1 %d  %d  \n",num_dims1, dims1[i]);
    }
    */

    float max_dist,fac,dist,dx,dy,bna,bua,descd,ov;
    for(int j=0,f1=0;j<dims1[1];j++,f1+=dims1[0]){    
      max_dist=sqrt(feat1[f1+5]*feat1[f1+6]);
      if(common_part)fac=30/max_dist;
      else fac=3;
      max_dist=max_dist*4;
      fac=1.0/(fac*fac);
      feat1a[2]=fac*feat1[f1+2];
      feat1a[3]=fac*feat1[f1+3];
      feat1a[4]=fac*feat1[f1+4];
      feat1a[7] = sqrt(feat1a[4]/(feat1a[2]*feat1a[4] - feat1a[3]*feat1a[3]));
      feat1a[8] = sqrt(feat1a[2]/(feat1a[2]*feat1a[4] - feat1a[3]*feat1a[3]));
      for (int i=0,f2=0; i<dims2[1]; i++,f2+=dims1[0]){
        //compute shift error between ellipses
        dx=feat2[f2]-feat1[f1];
        dy=feat2[f2+1]-feat1[f1+1];
        dist=sqrt(dx*dx+dy*dy);
        if(dist<max_dist){
    feat2a[2]=fac*feat2[f2+2];
    feat2a[3]=fac*feat2[f2+3];
    feat2a[4]=fac*feat2[f2+4];
    feat2a[7] = sqrt(feat2a[4]/(feat2a[2]*feat2a[4] - feat2a[3]*feat2a[3]));
    feat2a[8] = sqrt(feat2a[2]/(feat2a[2]*feat2a[4] - feat2a[3]*feat2a[3]));
    //find the largest eigenvalue
    float maxx=ceil((feat1a[7]>(dx+feat2a[7]))?feat1a[7]:(dx+feat2a[7]));
    float minx=floor((-feat1a[7]<(dx-feat2a[7]))?(-feat1a[7]):(dx-feat2a[7]));
    float maxy=ceil((feat1a[8]>(dy+feat2a[8]))?feat1a[8]:(dy+feat2a[8]));
          float miny=floor((-feat1a[8]<(dy-feat2a[8]))?(-feat1a[8]):(dy-feat2a[8]));

  float mina=(maxx-minx)<(maxy-miny)?(maxx-minx):(maxy-miny);
  float dr=mina/50.0;
  bua=0;bna=0;int t1=0,t2=0;
  //compute the area
  for(float rx=minx;rx<=maxx;rx+=dr){
    float rx2=rx-dx;t1++;
    for(float ry=miny;ry<=maxy;ry+=dr){
      float ry2=ry-dy;
      //compute the distance from the ellipse center
      float a=feat1a[2]*rx*rx+2*feat1a[3]*rx*ry+feat1a[4]*ry*ry;
      float b=feat2a[2]*rx2*rx2+2*feat2a[3]*rx2*ry2+feat2a[4]*ry2*ry2;
      //compute the area
      if(a<1 && b<1)bna++;
      if(a<1 || b<1)bua++;
    }
  }
  ov=100.0*(1-bna/bua);
  tover_out[j*dims2[1]+i]=ov;
  mover_out[j*dims2[1]+i]=ov;
  //printf("overlap %f  \n",over_out[j*dims2[1]+i]);return;
      }else {
  tover_out[j*dims2[1]+i]=100.0;
  mover_out[j*dims2[1]+i]=100.0;
      }
      descd=0;
      for(int v=9;v<dims1[0];v++){
  descd+=((feat1[f1+v]-feat2[f2+v])*(feat1[f1+v]-feat2[f2+v]));
      }
      descd=sqrt(descd);
      tdesc_out[j*dims2[1]+i]=descd;  
      mdesc_out[j*dims2[1]+i]=descd;  
    }
  }

    float minr=100;
    int mini=0;
    int minj=0;
    do{
        minr=100;
        for(int j=0;j<dims1[1];j++){    
    for (int i=0; i<dims2[1]; i++){
      if(minr>tover_out[j*dims2[1]+i]){
        minr=tover_out[j*dims2[1]+i];
        mini=i;
        minj=j;
      }
    }
        }
        if(minr<100){
    for(int j=0;j<dims1[1];j++){
      tover_out[j*dims2[1]+mini]=100;
    }   
    for (int i=0; i<dims2[1]; i++){
      tover_out[minj*dims2[1]+i]=100;
    }
    over_out[minj*dims2[1]+mini]=minr;
        }
    }while(minr<70);

    int dnbr=0;
    do{
      minr=1000000;
      for(int j=0;j<dims1[1];j++){    
        for (int i=0; i<dims2[1]; i++){
    if(minr>tdesc_out[j*dims2[1]+i]){
      minr=tdesc_out[j*dims2[1]+i];
      mini=i;
      minj=j;
    }
        }
      }
      if(minr<1000000){
        for(int j=0;j<dims1[1];j++){
    tdesc_out[j*dims2[1]+mini]=1000000;
    }   
        for (int i=0; i<dims2[1]; i++){
    tdesc_out[minj*dims2[1]+i]=1000000;
        }
        desc_out[minj*dims2[1]+mini]=dnbr++;//minr
      }
    }while(minr<1000000);

    delete []odim;
    delete []tdesc_out;
    delete []tover_out;
    delete []feat1a;
    delete []feat2a;
  }

Here is the memory profiler output:

Many thanks

Matthew

[James Tursa]

This is too much code to sift through. Can you simplify it to a much smaller size that still reproduces the problem?

[Matthew Gillman]

Hi there.

Yes, sorry, I realise it was too much :-)

I refactored my code so I had a perl wrapper which continuously called my MATLAB code (for one image), then exited MATLAB so as to free up memory, then called it again, etc. But I still got a memory allocation error after a while.

I realised that the problem was not a memory leak (where memory is gradually lost over time). Rather, the problem was with the size of the matrices that MATLAB was trying to allocate. The size of these matrices depends on files which specify the number of "feature points" (points of interest) in an image. It couldn't handle a file with 7678 points, but managed 6552, so I am guessing the limit is somewhere between these two, and basically to get round it I will have to only process those files which do not have too many feature points in. This should stop the memory allocation errors.

Thanks

Matthew

Answer 2

[Jan]

Are you sure that the 3rd input of the Mex function is a double?

double const *flag = (double *)mxGetData(prhs[2]);

In the caller it is mentioned as a flag with the values 0 and 1. It could be a logical also. So prefer:

double flag = mxGetScalar(prhs[2]);

It is better to rely on mwSize instead of int:

int const *dims1 = mxGetDimensions(prhs[0]);

mxGetDimension does not reply int.

But I do not think that this solves your problem.