Delaunay command gives fewer of triangles than expected in MATLAB

Question

I have got two unregistered images and a base image I use as reference for registration, image registration is performed as demonstrated in matlab example using SURF, now I have all images of 100*100 so after applying transformation matrix on both and saving all coordinates from three images in matrix named registeredPts when I apply delaunay command on a 30,000x2 matrix I get only 20,000 triangles where as as far as I know i should get 60,000 triangles approx

I have to use delaunay triangulation for image interpolation. I cannot figure out why so few number of triangles are formed though i cannot find any fault in feature based registration

close all
clear all
K = 2;
P1 = imread('C:\Users\Javeria Farooq\Desktop\project images\a.pgm');
%apply a command here that rgb2gray if it is a if it is rgb so convert it
% P1=rgb2gray(P1);
%reads the image to be registered
P2 = imread('C:\Users\Javeria Farooq\Desktop\project images\b.pgm');
% P2=rgb2gray(P2);

P3 = imread('C:\Users\Javeria Farooq\Desktop\project images\c.pgm');
% P3=rgb2gray(P3);
%reads the base image
 image1_gray = makelr(P1, 1, 100, 1/2);
%image1_gray = P1;
% makes lr image of first
 image2_gray= makelr(P2, 1, 100, 1/2);
image3_gray= makelr(P3, 1, 100, 1/2);
%image2_gray= P2;
%makes lr image of second
figure(1),imshow(image1_gray)
axis on;
grid on;
title('Unregistered image');
figure(2),imshow(image3_gray)
axis on;
grid on;
title('Unregistered image2');
figure(3),imshow(image2_gray)
axis on;
grid on;
title('Base image ');
impixelinfo
% both image displayed with pixel info
hold on
points_image1= detectSURFFeatures(image1_gray, 'NumScaleLevels', 100, 'NumOctaves', 12,  'MetricThreshold', 500 );
%detects surf features of first image
points_image2 = detectSURFFeatures(image2_gray, 'NumScaleLevels', 100, 'NumOctaves', 12,  'MetricThreshold', 500 );
points_image3 = detectSURFFeatures(image3_gray, 'NumScaleLevels', 100, 'NumOctaves', 12,  'MetricThreshold', 500 );
%detects surf features of second image
[features_image1, validPoints_image1] = extractFeatures(image1_gray, points_image1);
[features_image2, validPoints_image2] = extractFeatures(image2_gray, points_image2);
[features_image3, validPoints_image3] = extractFeatures(image3_gray, points_image3);
%extracts features of both images 
indexPairs = matchFeatures(features_image1, features_image2, 'Prenormalized', true) ;
indexPairs1 = matchFeatures(features_image3, features_image2, 'Prenormalized', true) ;
 % get matching points
matched_pts1 = validPoints_image1(indexPairs(:, 1));
matched_pts2 = validPoints_image2(indexPairs(:, 2));
matched_pts3 = validPoints_image3(indexPairs1(:, 1));
matched_pts4=validPoints_image2(indexPairs1(:, 2));
figure(4); showMatchedFeatures(image1_gray,image2_gray,matched_pts1,matched_pts2,'montage');
legend('matched points 1','matched points 2'); 
figure(5); showMatchedFeatures(image3_gray,image2_gray,matched_pts3,matched_pts4,'montage');
%matched features of both  images are displayed
legend('matched points 3','matched points 2'); 
 % Compute the transformation matrix using RANSAC
[tform, inlierFramePoints, inlierPanoPoints, status] = estimateGeometricTransform(matched_pts1, matched_pts2, 'projective');
[tform1, inlierFramePoints, inlierPanoPoints, status] = estimateGeometricTransform(matched_pts3, matched_pts4, 'projective');
%figure(6); showMatchedFeatures(image1_gray,image2_gray,inlierPanoPoints,inlierFramePoints,'montage');
%tform = estimateGeometricTransform(matched_pts1,matched_pts2,'projective')
%calculate transformation matrix using projective transform
T=tform.T;
r=[];
A=[];
l=1
[N1 N2]=size(image2_gray)
registeredPts = zeros(N1*N2,2);
pixelVals = zeros(N1*N2,1);
for row = 1:N1
    for col = 1:N2
        pixNum = (row-1)*N2 + col;
        pixelVals(pixNum,1) = image2_gray(row,col);
        registeredPts(pixNum,:) = [col,row];
    end
end
[r]=transformPointsForward(tform,registeredPts(:,:));
[q]=transformPointsForward(tform1,registeredPts(:,:));
%coordinates of base image
image2_gray=double(image2_gray);
R=2;
r1=r(:,1);
r2=r(:,2);
for row = 1:N1
        for col = 1:N2
            pixNum = N1*N2 + (row-1)*N2 + col;
            pixelVals(pixNum,1) = image1_gray(row,col);
            registeredPts(pixNum,:) = [r1(col,:),r2(col,:)];
        end
end
q1=q(:,1);
q2=q(:,2);
for row = 1:N1
        for col = 1:N2
            pixNum = N1*N2 +N1*N2+ (row-1)*N2 + col;
            pixelVals(pixNum,1) = image3_gray(row,col);
            registeredPts(pixNum,:) = [q1(col,:),q2(col,:)];
        end
end





 tri = delaunayTriangulation();
 tri.Points=[registeredPts(:,1),registeredPts(:,2)];

Answer 1

i figured out the problem in code

for row = 1:N1
        for col = 1:N2
            pixNum = N1*N2 + (row-1)*N2 + col;
            pixelVals(pixNum,1) = image1_gray(row,col);
            registeredPts(pixNum,:) = [r1(col,:),r2(col,:)];
        end
end

here same registered points are reassigned again and again, i corrected it by assigning the points out of for loop using registeredPts=[registeredPts;r]

now i m getting above 60,000 triangles which works fine