classdef ImgProcess < TestImgGen & ImgProcess_CheckerBoard
    %IMGPROCESS Summary of this class goes here
    %   Detailed explanation goes here
    
    properties
        
        CamPlane
        CamPlaneCalibrationSaved
        
        darkfieldImg
        clearfieldImg
        darkFieldStruct
        clearFieldStruct
        
        calibOutData
        
        
        
        
    end
    
    methods
        function obj = ImgProcess()
            %IMGPROCESS Construct an instance of this class
            %   Detailed explanation goes here
            obj.CamPlane = [];
            obj.CamPlaneCalibrationSaved = false;
            obj.testPatterns([360,640,3],[20 20],2,1);
        end
        
        function borderOut = detectBorder(obj,img)
            img = rescale(img);
            binaryBorder = imfill(imbinarize(img,'global'),'holes');
            borderOut = bwboundaries(binaryBorder);
        end
        function [centroid,scale] = fitBiggestRectangle(obj,im_size,border_points)
            xy = regionprops(border_points,'PixelList');
            xy = cat(1, xy.PixelList);
            centroid = mean(xy);
            
            rectangle = flip(im_size(1:2));
            L = 0;
            U = 1;
            
            while L/U < 0.99998
                scale = (L+U)/2;
                if obj.iscolision(centroid, rectangle, scale, xy)
                    U = scale;
                else
                    L = scale;
                end
                
            end
            %function point_in_rectangle()
        end
        function ret = iscolision(~,center, side, scale, xy)
            x1 = center(2) - side(2)*scale/2;
            x2 = center(2) + side(2)*scale/2;
            
            y1 = center(1) - side(1)*scale/2;
            y2 = center(1) + side(1)*scale/2;
            
            x = xy(:,2);
            y = xy(:,1);
            if any((x > x1) & (x < x2) & (y > y1) & (y < y2))
                ret = true;
            else
                ret = false;
            end
        end
        
        
%         function [px,py,pz]=projection(~,a,b,c,d,x,y,z)
%             % given an plane equation ax+by+cz=d, project points xyz onto the plane
%             % return the coordinates of the new projected points
%             % written by Neo Jing Ci, 11/7/18
%             A=[1 0 0 -a; 0 1 0 -b; 0 0 1 -c; a b c 0];
%             B=[x; y; z; d];
%             X=A\B;
%             px=X(1);
%             py=X(2);
%             pz=X(3);
%         end
%         
%         function [X,Y,Z] = rectangle4(~,x,y,h,w)
%             h = h-1; w = w-1;
%             X = [x;x;x+w;x+w];
%             Y = [y;y+h;y+h;y];
%             Z = [0;0;0;0];
%         end
%         
%         function [imw_crop, imw_same] = changePerspective(~,img,yaw,pitch,roll,camPos,F)
%             % Source:
%             % https://www.mathworks.com/matlabcentral/answers/166497-create-a-projective-homography-matrix-with-pitch-roll-yaw
%             
%             % Set camera position here. To understand the effect of angles, move it
%             % along the z-axis and y-axis so it is above the centre of the top edge of
%             % the picture. Note that as cammoves is true in call to homography_matrix,
%             % the z-axis movement needs to be negative.
%             % The camera is looking along the z-axis.
%             % campos = [size(img,2)/2; 0; -200];
%             
%             % F = 100;    % focal length
%             
%             
%             yawmat = rotation_matrix([0 0 1], -yaw);       % about z
%             pitchmat = rotation_matrix([0 1 0], -pitch);   % about y
%             rollmat = rotation_matrix([1 0 0], -roll);     % about x
%             
%             rmat = rollmat * pitchmat * yawmat;   % zyx order - default for angle2dcm
%             
%             % construct homography matrix
%             hmat = homography_matrix([], rmat, camPos, F, true);
%             
%             % construct projective2d object
%             tform = projective2d(hmat.');
%             
%             %pad image
%             %             sImg = size(img);
%             %             if sImg(1) > sImg(2)
%             %                 img = [zeros(sImg(1),floor((sImg(1)-sImg(2))/2),sImg(3)),...
%             %                     img,...
%             %                     zeros(sImg(1),floor((sImg(1)-sImg(2))/2),sImg(3))];
%             %             else
%             %                img = [zeros(ceil((sImg(2)-sImg(1))/2),sImg(2),sImg(3));...
%             %                     img;...
%             %                     zeros(ceil((sImg(2)-sImg(1))/2),sImg(2),sImg(3))];
%             %             end
%             
%             
%             % image warped and cropped
%             imw_crop = imwarp(img, tform);
%             
%             % image warped in same coordinate system as original
%             imw_same = imwarp_same(img, tform);
%         end
        
       
        
        %% QR Position Marker detection
        
        function [bwHollow, bwEnhanced, bwBinarizedOnly, RGBvisual] = AdaptiveBinarizeAndFilterHollowObjects(~,img,ImageBinarizationSensitivity,QRMarkLineThickness)
            % AdaptiveBinarizeFilterHollowObjects Binarize image (RGB or
            % grayscale) with adaptive threshold. Filter Hollow objects
            
            % RGBvisual  ==1 (bright) bwHollow
            %            ==0.5 (dimmed) ~bwHollow
            %            green = added by enhancement
            %            red = removed by enhancement
            
            % image input normalization
            if size(img,3) > 1
                img = rgb2gray(img);
            end
            if ~isa(img,'double')
                img = im2double(img);
            end
            
            % estimate of average background illumination (sensitivity best resoults = 0)
            T = adaptthresh(img,ImageBinarizationSensitivity,'ForegroundPolarity','dark','Statistic','gaussian');
            
            % binarize image with adaptive treshold
            bwBinarizedOnly = ~imbinarize(img,T);
            
            % inverse
            bwo = bwBinarizedOnly;
            bwEnhanced = imclearborder(bwo);
            
            % clean binary image (reduces amount of objects to test)
            if QRMarkLineThickness > 0
                % supress objects with smaller thickness than (3)
                bwEnhanced = imopen(bwo, ones(QRMarkLineThickness));
                
                % fill small holes inside QR marks
                bwEnhanced = imclose(bwEnhanced, ones(QRMarkLineThickness));
            end
            
            % select binary shapes with hole(s) in them
            % Euler number of Binary image:
            % total number of objects in the image (1) minus the number of holes in those objects
            bwHollow = bwpropfilt(bwEnhanced, 'EulerNumber', [-Inf, 0]);
            
            
            % Binarization visualization (could be minimized)
            opacity = 0.5;
            
            pxRem = bwBinarizedOnly & ~bwEnhanced;
            pxAdd = bwEnhanced & ~bwBinarizedOnly;
            
            W = bwHollow & ~(pxRem | pxAdd);
            Gr = ~bwHollow & (bwBinarizedOnly | bwEnhanced) & ~(pxRem | pxAdd);
            
            R = (pxRem & bwHollow | W)*1 + (pxRem & ~bwHollow | Gr)*opacity ;
            
            G = (pxAdd & bwHollow | W)*1 + (pxAdd & ~bwHollow | Gr)*opacity;
            
            B = W*1 + Gr*opacity;
            
            
            RGBvisual = cat(3,R, G, B);
            
        end
        
        function [labeledMask, markersPosition] = findMarkers(obj,bwHollow,bwSource,QRMarkPatternPolicy,QRMarkPatternTolerance)
            
            bw2 = bwSource;
            bw3 = bwHollow;
            
            
            
            % binary mask filled hole in middle
            bmask = imfill(bw3,'holes');
            
            % get centre of those holes in binary shapes
            stats = regionprops(bmask & ~bw3,'Centroid');
            
            % create matrix labling binary components
            bwl = zeros([size(bmask),length(stats)]);
            
            flag = zeros(length(stats),1);
            for i = 1:length(stats)
                % position of centroid
                cx = round(stats(i).Centroid(1));
                cy = round(stats(i).Centroid(2));
                
                % find roi
                bwl(:,:,i) = bwselect(bmask,cx,cy);
                roi = bwl(:,:,i)&bw2;
                
              %  imshow(roi)
                flag(i) = obj.checkXY([cx,cy],roi,QRMarkPatternPolicy,QRMarkPatternTolerance);
                
              %  flag(i)
            end
            
            % output
            cID = find(flag==1);
            
            Centroids = cat(1, stats.Centroid);
            markersPosition = Centroids(cID,:);
            
            labeledMask = zeros(size(bwSource));
            
            if ~isempty(cID)
                % clockwise order
                K = convhull(markersPosition(:,1),markersPosition(:,2));
                K = K(1:end-1);
                markersPosition = markersPosition(K,:);
                [markersPosition(:,1),markersPosition(:,2)] = poly2cw(markersPosition(:,1),markersPosition(:,2));
                % first element closest to origin
                [~,K] = min(sum(markersPosition.^2, 2));
                markersPosition = circshift(markersPosition, (length(cID)-K+1), 1);
                
                % corresponding mask
                for i = 1:length(cID)
                    [r,~] = find( markersPosition(i,:) == Centroids );
                    labeledMask = labeledMask + i*(bwl(:,:,r(1))&bw2);
                end
                
            end
            
            
            
        end
        
        function isMarker = checkXY(obj,cp, roi, policy, tolerance)
            % center [x,y] coordinate, policy in format (1of1, 1of2, 2of2, 1of3, 2of3, 3of3)
            if length(policy)~=4
                error('Wrong "policy" input')
            end
            
            n = str2double(policy(1));
            m = str2double(policy(4));
            if isempty(n) || isempty(m) || n>m || n<1 || m>3
                error('Wrong "policy" input')
            end
            
            % diagonal endpoints y=kx+q, 45°->k=1;
            q = (cp(2) - cp(1));
            maxX = size(roi,2);
            if maxX + q > size(roi,1), maxX = size(roi,1) -q;end % fixed?
            d = improfile(roi,[-q, maxX],[0,maxX+q]);
            d(isnan(d))=[];
            d = logical(d);
            if isempty(d)
                error('Empty diagonal (should not happen)')
            end
            
            % vertical, horizontal, diagonal slice
            stripes = {roi(:,cp(1)), roi(cp(2),:)', d};
            
            isMarker = zeros(1,3,'logical');
            
            for i = 1:m
                if any(isnan(stripes{i})) || isempty(stripes{i})
                    error('');
                end
                
                isMarker(i) = obj.QRcompareBWRatio(stripes{i},tolerance);
                if sum(isMarker) >= n
                    isMarker = true;
                    return;
                end
            end
            
            isMarker = false;
        end
        
        function isMarker = QRcompareBWRatio(~,strip,tolerance)
            % modification of:
            % Kei Otsuka (2020). QR Code reader example 
            % (https://www.mathworks.com/matlabcentral/fileexchange/71090-qr-code-reader-example),
            % MATLAB Central File Exchange.
            
            % QR Marker (inverted) ...._?_???_?_.... , pulse width ratio [1 1 3 1 1]
            % Marker line width tolerance according to QR specifiaction = 0.5
            
            pw = pulsewidth(double(strip)); % width of pulses (==1)
            ps = pulsesep(double(strip)); % width of 0 between pulses
            
            isMarker = false;
            
            
            % if 3+ pulses with 2+ gaps and puseses == gaps + 1
            flagout = (size(ps, 1) >= 2 && size(pw, 1) >= 3 && (size(ps, 1)+1 == size(pw, 1)) );
            
           
            if flagout == 1
                for i = 0:length(pw)-3
                    pat = [pw(1+i), ps(1+i), pw(2+i), ps(2+i), pw(3+i)];
                    lWidth = median([pw(1+i), ps(1+i), pw(2+i)/3, ps(2+i), pw(3+i)]);
                    pat = pat/lWidth;
                    % compare thickness sequence with specification :
                    % (+-0.5) tolerance, allow one edge to be thicker
                    % (shadow)
                    if pdist2(pat(1:4), [1 1 3 1],'chebychev') < tolerance || pdist2(pat(2:end), [1 3 1 1],'chebychev') < tolerance
                        isMarker = true;
                    end
                end
            end
            
%             
%             
%             if flagout == 1
%                 width = median([pw(1), ps(1), pw(2)/3, ps(2), pw(3)]);
%                 pat = [pw(1), ps(1), pw(2), ps(2), pw(3)] / width;
%                 % compare thhickness sequence with specification : (+-0.5) tolerance
%                 if pdist2(pat(1:3), [1 1 3]) < tolerance^2*5 || pdist2(pat(3:end), [3 1 1]) < tolerance^2*5
%                     isMarker = true;
%                 end
%             end
        end
        
        function  [tform, refP, PlaneRectangle,imScale] = fitGeoTransformationPreservingPhysicalDim(obj,markersPosition)
            [refP,imScale, PlaneRectangle] = obj.GetRefferencePoints(markersPosition);
            tform = fitgeotrans(markersPosition,refP,'projective');
        end
        
        function [Pout,imScale, PlaneRectangle] = GetRefferencePoints(~,Pin)
            % imScale [px/mm]
            x1 = mean(Pin(1:2,1));
            x2 = mean(Pin(3:4,1));
            y1 = mean(Pin([1,4],2));
            y2 = mean(Pin(2:3,2));
            
            % physical distance between marks (mm)
            markerWidth = 15;
            xReal = 215.92-markerWidth;
            yReal = 166-markerWidth;
            
            % image scale
            xScale = (x2-x1)/xReal;
            yScale = (y2-y1)/yReal;
            
            TargetScale = mean([xScale,yScale]);
            
            % correct x1, x2 position
            xMove = (xReal*TargetScale-(x2-x1))/2;
            
            x1 = x1 - xMove;
            x2 = x2 + xMove;
            
            % correct y1, y2 position
            yMove = (yReal*TargetScale-(y2-y1))/2;
            
            y1 = y1 - yMove;
            y2 = y2 + yMove;
            
            Pout = [x1 y1;...
                x1 y2;...
                x2 y2;...
                x2 y1];
            
            imScale = TargetScale;
            PlaneRectangle = [x1-markerWidth/2*imScale, y1-markerWidth/2*imScale,...
                x2-x1+markerWidth*imScale, y2-y1+markerWidth*imScale];
        end
        
        function [outIm2, markersOutIm2, outIm1, markersOutIm1, rectOutIm1] = WarpAndCropImg(obj,imIn, Tform, markersP, cropRectangle)
            % Warp and crop image, transform detected points
            
            % use saved parameters if inly imIn
            if nargin == 2 
                if ~obj.CamPlaneCalibrationSaved
                    error("No Camera-Plane calibration saved");
                end
                Tform = obj.CamPlane.TForm;
                markersP = obj.CamPlane.QRMarkersPos;
                cropRectangle = obj.CamPlane.Border;               
            elseif nargin < 5
                    error("Not enough arguments");
            end
            
            % warp image
            [outIm1,RB] = imwarp(imIn,Tform);
            
            % transform points coordinates
            [x,y] = transformPointsForward(Tform,markersP(:,1),markersP(:,2));
            [xx,yy] = worldToIntrinsic(RB,x,y);
            markersOutIm1 = [xx,yy];
            
            % transform rectangle coordinates
            [rx,ry] = worldToIntrinsic(RB,[cropRectangle(1),cropRectangle(1)+cropRectangle(3)],[cropRectangle(2),cropRectangle(2)+cropRectangle(4)]);
            rectOutIm1 = [rx(1) ry(1), rx(2)-rx(1),ry(2)-ry(1)];
            
            % crop warped image
            outIm2 = imcrop(outIm1,rectOutIm1);
            xxx = xx - rectOutIm1(1);
            yyy = yy - rectOutIm1(2);
            markersOutIm2 = [xxx,yyy];            
        end
        
        %% brightness
        function [imgOut] = addCornersToImg(~,imgIn,sizeCorner,positionM)
            % create mask with 1 in top left
            mask = zeros([size(imgIn,1),size(imgIn,2)]);
            mask(1:sizeCorner(1),1:sizeCorner(2)) = 1;
            
            imgOut = imgIn;
            for i = 1:4
                i_mask = mask;
                % if bottom flip mask
                if mod(i,2)
                    i_mask = flip(i_mask,1);
                end
                % if right flip mask
                if i > 2
                    i_mask = flip(i_mask,2);
                end
                % eg. [1 0; 0 -1]
                if positionM(i) == 1
                    % white corner
                    imgOut = imgOut .* ~i_mask + i_mask*1;
                elseif positionM(i) == -1
                    % black corner
                    imgOut = imgOut .* ~i_mask + i_mask*0;
                end
            end           
        end
        
        function [imStruct,imCombined,HistVal,mask, c1, c2, c3, c4] = findCornersInImg(obj,imStruct, imMasksCorner, c1, c2, c3, c4)
            % find projected ref points to normalize snapshot of projected image
            % img expected to be WarpedAndCropedImg 
            G = 10;
            for i = 1:length(imStruct)
                imStruct(i).imgGray = im2double(rgb2gray(imStruct(i).img));
            end
            if nargin < 3
             [imStruct(1).mask,c1, imStruct(2).mask,c2] = obj.findBrightSpot(imStruct(1).imgGray,...
                 imStruct(2).imgGray,G);
             [imStruct(3).mask, c3, imStruct(4).mask, c4] = obj.findBrightSpot(imStruct(3).imgGray,...
                 imStruct(4).imgGray,G);
            end     
             % save values of dark & bright spots to struct
             bright = [ median(imStruct(1).imgGray(imStruct(1).mask)),...
                 median(imStruct(2).imgGray(imStruct(2).mask)),...
                 median(imStruct(3).imgGray(imStruct(3).mask)),...
                 median(imStruct(4).imgGray(imStruct(4).mask))];
             
             dark = [ median(imStruct(1).imgGray(imStruct(2).mask)),...
                 median(imStruct(2).imgGray(imStruct(1).mask)),...
                 median(imStruct(3).imgGray(imStruct(4).mask)),...
                 median(imStruct(4).imgGray(imStruct(3).mask))];
         
             % centroid -> quadrant to mask
             sImg = size(imStruct(1).imgGray);
             q12 = [c1(2) < sImg(1)/2, c1(2) > sImg(1)/2] & [c1(1) < sImg(2)/2; c1(1) > sImg(2)/2] |...
                 [c2(2) < sImg(1)/2, c2(2) > sImg(1)/2] & [c2(1) < sImg(2)/2; c2(1) > sImg(2)/2];
             
             q34 = [c3(2) < sImg(1)/2, c3(2) > sImg(1)/2] & [c3(1) < sImg(2)/2; c3(1) > sImg(2)/2] |...
                 [c4(2) < sImg(1)/2, c4(2) > sImg(1)/2] & [c4(1) < sImg(2)/2; c4(1) > sImg(2)/2];
             % combine imgs binary mask find
             Q12 = imresize(q12,[sImg],'nearest');
             Q34 = imresize(q34,[sImg],'nearest');
             % temporary images
             tImg1 = imStruct(1).imgGray;
             tImg2 = imStruct(3).imgGray;
             % mask images with NaN
             tImg1(Q12) = NaN;
             tImg2(Q34) = NaN;
             % temporary img 
             tImg(:,:,1) = tImg1;
             tImg(:,:,2) = tImg2;
             % mean tImg (x,y,2) -> (x,y,1)
             imCombined = nanmean(tImg,3) ;
             HistVal = [bright;dark];
             mask = (imStruct(1).mask | imStruct(2).mask | imStruct(3).mask | imStruct(4).mask);
        end
        
        function [imMask1, c1, imMask2, c2] = findBrightSpot(~,im1,im2,G)
            i1 = im1-im2; i1(i1<0) = 0;
            i1 = imadjust(i1);
            
            i2 = im2 - im1; i2(i2<0) = 0;
            i2 = imadjust(i2);
            
            T1 = adaptthresh(imcomplement(i1),0.5,'ForegroundPolarity','dark');
            T2 = adaptthresh(imcomplement(i2),0.5,'ForegroundPolarity','dark');
            
            imMask1 = bwareafilt(imbinarize(imadjust(i1),T1),1);
            imMask2 = bwareafilt(imbinarize(imadjust(i2),T2),1);
            
            stats = regionprops(imMask1,'Centroid');
            c1 = stats.Centroid;
            
            stats = regionprops(imMask2,'Centroid');
            c2 = stats.Centroid;
            
        end
        
        %% checker board - majority moved to ImgProcess_CheckerBoard subclass    
        function [I, stdValBg] = AdjustAndCropCapturedImg(obj,I,bcg,allowImAdjust,litAreaBorderThickness)
            imOut = I;
            
            % dark field & clearfield load
            dd =   obj.darkfieldImg;
            stdValBg = std(dd(:));
            cc =  obj.clearfieldImg;
            
            % estimante white and black levels, match to clearfield &
            % darkfield levels
            whMask = imbinarize(imOut-dd);
            litMask = imfill(imbinarize(cc-dd),'holes');
            
            wL = median(imOut(whMask & litMask));
            dL = median(imOut(~whMask & litMask));
            
            cfL = median(cc(:));
            dfL = median(dd(:));
            
            m = (cfL-dfL)/(wL-dL);
            
            % adjust histogram of captured img
            imOut = imOut*m - ( dL*m -  dfL );
            %imOut(imOut<0) = 0; imOut(imOut>1) = 1;
            
            % adaptive brigthness correction
            imOut2 =  ((imOut-dd)./(cc-dd)) .* litMask;
            
            w = max(imOut2(:)) - min(imOut2(:));                
            imOut2 = imOut2./w - (min(imOut2(:))./w);
            if allowImAdjust
                imOut2 = imadjust(imOut2);
            end
            imOut2(~litMask) = NaN;
            I = imOut2; 
            
            % final image adjust - crop nan rows, remaining nan turn to
            % whine/back
            I(I>1) = 1;
            nanMask = isnan(I);
            if nargin <5, litAreaBorderThickness = 3;end
            nanMask = imdilate(nanMask,strel('sphere',round(litAreaBorderThickness/640*length(imOut))));
            I(nanMask) = NaN;
            
            if nargin> 2 && bcg == 1
                NaNPxWalue = median(I(whMask & litMask),'omitnan');
            else
                NaNPxWalue = median(I(~whMask & litMask),'omitnan');
            end
            
            % create nan crop mask
            if isempty(obj.calibOutData) || isempty(obj.calibOutData.NaNCropmask)
                obj.calibOutData.NaNCropmask = ~any(~isnan(I), 2);
            end
          
            I = obj.cropNaNAndResize(I);
            I(isnan(I)) = NaNPxWalue; 
        end
        
        function I = cropNaNAndResize(obj,I)
            if isempty(obj.calibOutData) || isempty(obj.calibOutData.NaNCropmask)
                error('Run AdjustAndCropCapturedImg first')
            end
              I(obj.calibOutData.NaNCropmask,:)=[];   
            % resize
            %I = imresize(I,640/length(I),'Method','nearest');
             if any(size(I) ~= [360 640])
                imgScales = [360 640]./size(I);
                imgScale = min(imgScales(1:2));
                I = imresize(I, round(imgScale*size(I)),'Method','nearest');
             end
            
             IR=I;
                imgSizeR = size( IR);
                toFill = [360 640] - imgSizeR;
                I = ones(360, 640).*NaN;
                I(1+floor(toFill(1)/2):360-ceil(toFill(1)/2),...
                    1+floor(toFill(2)/2):640-ceil(toFill(2)/2)) =  IR(:,:);
            
        end
        
    end
end


