title: super resolution of low resolution images using delaunay triangulation, negative pixel values for the resultant High resolution image
body: <p>i have to do super resolution of two low resolution images to obtain a high resolution image.</p>
  
  <p>2nd image is taken as base image and the first image is registered with respect to it . i used SURF algorithm for image registration . A Delaunay triangulation is
  constructed over the points using a built-in MATLAB delaunay
  function . The HR grid of size is
  constructed for a prespecified resolution enhancement factor R Then HR algorithm for interpolating the pixel values on the
  HR grid is summarized next.</p>
  
  <p>HR Algorithm Steps:
  1. Construct the Delaunay triangulation
  over the set of scattered vertices in the
  irregularly sampled raster formed from the
  LR frames.</p>
  
  <ol>
  <li><p>Estimate the gradient vector at each
  vertex of the triangulation by calculating the unit normal vector of neighbouring vector using cross product method.Sum of the unit normal vector of each triangle  multiplied by its area is divided by summation of area of all neighbouring triangles to get the vertex normal.</p></li>
  <li><p>Approximate each triangle patch in
  the triangulation by a continuous and,
  possibly, a continuously differentiable
  surface, subject to some smoothness constraint.
  Bivariate polynomials or splines
  could be the approximants as explained
  below.</p></li>
  <li>Set the resolution enhancement factor
  along the horizontal and vertical directions
  and then calculate the pixel value
  at each regularly spaced HR grid point to
  construct the initial HR image </li>
  </ol>
  
  <p>The bivariate polynomial i used is mentioned in the code, using pixel values at each vertex of a triangle and corresponding gradient in x and y directions i calculated the nine constants associated with each triangle then defined a high resolution grid , calculated the pixel values at each point using the constants calculated</p>
  
  <p>i am attaching my code with it, the problem i am facing is that i am just getting a gray image as out put HR image , because the constants i have calculated have negative values resulting in negative pixel values</p>
  
  <p>another problem i realized with my code is in gradient estimation i get a lot of 'NaN' as a result of gradient calculation. </p>
  
  <p>if any one can please spent some time to help me out</p>
  
  <pre><code>close all
  clear all
  K = 2;
  P1 = imread('C:\Users\Javeria Farooq\Desktop\project images\a.pgm');
  %reads the image to be registered
  P2 = imread('C:\Users\Javeria Farooq\Desktop\project images\b.pgm');
  %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);
  %image2_gray= P2;
  %makes lr image of second
  figure(1),imshow(image1_gray)
  axis on;
  grid on;
  title('Unregistered image');
  figure(2),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 );
  %detects surf features of second image
  [features_image1, validPoints_image1] = extractFeatures(image1_gray, points_image1);
  [features_image2, validPoints_image2] = extractFeatures(image2_gray, points_image2);
  %extracts features of both images 
  indexPairs = matchFeatures(features_image1, features_image2, 'Prenormalized', true) ;
   % get matching points
  matched_pts1 = validPoints_image1(indexPairs(:, 1));
  matched_pts2 = validPoints_image2(indexPairs(:, 2));
  figure; showMatchedFeatures(image1_gray,image2_gray,matched_pts1,matched_pts2,'montage');
  %matched features of both  images are displayed
  legend('matched points 1','matched points 2'); 
   % Compute the transformation matrix 
  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);
  % s= zeros(N1*N2,2);
  pixelVals = zeros(N1*N2,1);
  [N1 N2]=size(image2_gray)
  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);
  %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(row,1),r2(row,1)];
          end
      end
  % all pixel values are saved in pixelVals
  %all registered points are saved first base image then unregistered image
  %delaunay triangulation of all coordinates passing x and y coordinates from registered Points 
   tri = delaunay(registeredPts(:,1),registeredPts(:,2));
  figure(3), triplot(tri,registeredPts(:,1),registeredPts(:,2))
  save tri
  % Estimate the gradient vector at each vertex
  [totalTris,three] = size(tri);
  [totalPoints,two] = size(registeredPts);
  vGradientVecs = zeros(totalPoints,2);
  triAreas = zeros(totalTris,1);
  triUnitNormals = zeros(totalTris,3);
  vUnitNormals = zeros(totalPoints,3);
  %   1. Find the unit normal vectors and the areas of all triangles,
  %      then find the product of these two numbers for each triangle
  for triNum = 1:totalTris
      v = tri(triNum,:);
      % 3D triangle points: x,y,pixel
      b=pixelVals(v);
      b=b(:);
      p = [registeredPts(v,:),b];
      % triangle area
      triAreas(triNum) = polyarea([p(:,1)],[p(:,2)]);
  
      % directional vectors representing the surface of the plane
      d1 = p(2,:)-p(1,:);
      d2 = p(3,:)-p(1,:);
  
      % cross product of these vectors
      crossp = cross(d1,d2);
  
  % If u = [u1 u2 u3] and v = [v1 v2 v3], we know that the product w is defined as w = [(u2v3 – u3v2) (u3v1 - u1v3) (u1v2 - u2v1)]
      % normalized cross product = unit normal vector for the triangle
      dist = sqrt(sum(crossp.^2));
      triUnitNormals(triNum,:) = crossp./dist;
  end
  
  %   %2. %Estimate the unit normal vector at each vertex
  %       a. Find the triangle patches that neighbor the vertex
  %       b. Find the unit normal vectors of these regions
  %       c. Multiply each of these vectors by the area of the 
  %          associated region, then sum these numbers and divide
  %          by the total area of all the regions
  
   for pointNum = 1:totalPoints
      [neighbors,x] = find(tri==pointNum);
      areas = triAreas(neighbors);
      areas3 = [areas,areas,areas];
      triNormsSum = sum(triUnitNormals(neighbors,:).*areas3);
      triAreasSum = sum(areas);
       vUnormalized = triNormsSum./triAreasSum;
          vUnitNormals(pointNum,:) = ...
          vUnormalized./sqrt(sum(vUnormalized.^2));
      if( triAreasSum == 0 )
          triAreasSum = 0.0001;
      vUnormalized = triNormsSum./triAreasSum;
  
      % re-normalize
      vUnitNormals(pointNum,:) = ...
          vUnormalized./sqrt(sum(vUnormalized.^2));
  end
  
  %   3. Find the gradients along the x and y directions for each vertex
  %      vertex's unit normal: n = [nx,ny,nz]
  %      x-direction gradient: dz/dx = -nx/nz
  %      y-direction gradient: dz/dy = -ny/nz
  % 
   for pointNum = 1:totalPoints
      nz = vUnitNormals(pointNum,3);
      if( nz == 0 )
          nz = 0.0001;
      end
      vGradientVecs(pointNum,1) = -vUnitNormals(pointNum,1)./nz;
      vGradientVecs(pointNum,2) = -vUnitNormals(pointNum,2)./nz;
  % end
   end
   end
  %   1. Find the 3 equations for each vertex, and
  %      place them in c_equations matrix;
  % c_equations = [A for vertex 1;
  %                A for vertex 2; ...
  %                A for vertex totalPoints]
  % c(point,row,:) gives one row from an A matrix
   Btotal = zeros(3,totalPoints);
   c_equations = zeros(3*totalPoints,3,9);
   for pointNum = 1:totalPoints
  %     % B = [pixVal; x gradient; y gradient] at this vertex
   z = pixelVals(pointNum);
       B = [z; vGradientVecs(pointNum,1); vGradientVecs(pointNum,2)];
  %     
  %     % Compile all B matrices into a vector
       Btotal(:,pointNum) = B;
  
      % B = Ac to calculate c which is c=[c1 c2 .....c9]' take invA and
      % multiply by B
      x = registeredPts(pointNum,1);
      y = registeredPts(pointNum,2);
      A = [1   x   y   x^2  y^2  x^3     (x^2)*y  x*(y^2)  y^3; ...
           0   1   0   2*x   0   3*(x^2)  2*x*y   y^2      0; ...
           0   0   1   0    2*y  0        x^2     2*x*y    3*(y^2)];
  
      % Compile all A matrices into a vector
      c_equations(pointNum,1,:) = A(1,:);
      c_equations(pointNum,2,:) = A(2,:);
      c_equations(pointNum,3,:) = A(3,:);
  end
  
  %   2. Find the c values for each triangle patch
  c = zeros(totalTris,9);
  c9 = zeros(9,9);
  for triNum = 1:totalTris
      p1 = tri(triNum,1);
      p2 = tri(triNum,2);
      p3 = tri(triNum,3);
  
      B9 = [Btotal(:,p1); Btotal(:,p2); Btotal(:,p3)];
      c9 = [(c_equations(p1,1,:)); (c_equations(p1,2,:)); (c_equations(p1,3,:)); ...
            (c_equations(p2,1,:)); (c_equations(p2,2,:));( c_equations(p2,3,:)); ...
            (c_equations(p3,1,:)); (c_equations(p3,2,:));( c_equations(p3,3,:))];
        C9=squeeze(c9);
      c(triNum,:) = pinv(C9)*B9;   %linsolve(c9,B9);
  
  end
  % xc = findBPolyCoefficients1(tri,registeredPts,pixelVals,vGradientVecs);
  % save xc
  % %   2. For each point on the HR grid, find the associated triangle patch,
  % %      extract its c values, and use these values as the coefficients
  % %      in a bivariate polynomial to calculate the HR pixel value at
  % %      each grid point (x,y)
  [N1,N2]=size(image1_gray);
  [totalTris,three] = size(tri);
  
  M = N1*R-1;
  N = N2*R-1;
  HRimage = zeros(M,N);
  HRtriangles = zeros(M,N);
  [X,Y] = meshgrid(1:1/R:N2,1:1/R:N1);
  
  % Check all the triangles in order noting in which triangle each HR
  % grid point occurs.
  for triNum = 1:totalTris
      pts = registeredPts(tri(triNum,:),:);
      IN = inpolygon(X,Y,pts(:,1),pts(:,2));  % NxM
      HRtriangles(ind2sub(size(IN),find(IN==1))) = triNum;
  end
  % there is  a problem with this part of code , 
  for y = 1:M % row
      for x = 1:N   % col  
          % For testing, average the pixels from the vertices of the
          % triangle the HR point is in.
  %         pix = pixelVals(tri(HRtriangles(x,y),:));
     %      HRimage(x,y) = (pix(1) + pix(2) + pix(3))/3;
          % Extract appropriate set of 9 c values
           HRptC = c(HRtriangles(x,y),:);
           % Bivariate polynomial
      HRimage(x,y) = sum(HRptC.*[1,x,y,x^2,y^2,x^3,(x^2)*y,x*(y^2),y^3]);
           g(x,y)=HRimage(x,y);
          %changd xy with yx
      end
  end
  %   HRimage = estimateGridVals1(tri,registeredPts,R,N1,N2,pixelVals);
  % %Estimating Grid values at each patch
  %  %save HRimage
  
       g(g(:,:)&lt;0)=0;
  figure(8),imshow(g,[]);
  </code></pre>