exudates

exudates.cpp

#include <iostream>
#include <fstream>
#include "opencv2/opencv.hpp"
#include "opencv2/ximgproc.hpp"
using namespace cv;
using namespace cv::ximgproc;
using std::cout;
using std::endl;

// https://ieeexplore.ieee.org/document/8015123


// https://en.wikipedia.org/wiki/YIQ#Formulas
// this expects rgb order !
void enhance(const Mat &in, Mat &out) {
    Mat_<float> to_yiq(3,3);
    to_yiq << 0.299, 0.587, 0.144,
           0.596, -0.274, -0.322,
           0.211, -0.523, 0.312;
    transform(in, out, to_yiq);

    Mat chan[3];
    split(out, chan);
    Ptr<CLAHE> cl = createCLAHE(10);
    chan[0].convertTo(chan[0], CV_8U, 255);
    cl->apply(chan[0], chan[0]);
    chan[0].convertTo(chan[0], CV_32F, 1.0/255);
    merge(chan, 3, out);

    Mat_<float> from_yiq(3,3);
    from_yiq << 1, 0.956, 0.619,
           1, -0.272, -0.647,
           1, -1.106, 1.703;
    transform(out, out, from_yiq);
}



struct cluster {
    Mat features;
    Point center;
    int label;
};


// https://ieeexplore.ieee.org/document/8015123#Figure5
double contextual(std::vector<cluster> &clusters, int centerindex, double radius, double meangray) {
    Point p = clusters[centerindex].center;
    int num = 0;
    double ctxvalue = 0;
    double val = clusters[centerindex].features.at<double>(0); // mean_R
    for (size_t i=0; i<clusters.size(); i++) {
        if (i == centerindex)
            continue;
        Point p2 = clusters[i].center;
        double distsqr = sqrt((p.x-p2.x) * (p.x-p2.x)  + (p.y-p2.y) * (p.y-p2.y));
        if (distsqr < radius) {
            double pv = clusters[i].features.at<double>(0); // mean_R
            ctxvalue += sqrt((pv-val)*(pv-val)) / distsqr;
            num ++;
        }
    }
    return (val * ctxvalue) / (num * meangray);
}

int make_dataset()
{
    String diaret = "C:\\data\\diaret\\";
    std::vector<String> fn;
    glob(diaret + "images\\", fn);
    cout << fn.size() << endl;

    int min_element_size = 30;
    int region_size = 30;
    int num_iterations = 3;
    int algorithm = SLICO;
    float ruler = 0.1;

    FileStorage fs("exudates.xml", 1);
    fs << "exudates" << "[";

    for (size_t i=0; i<fn.size(); i++) {
        fs << "{:";
        fs << "image" << fn[i];
        Mat img = imread(fn[i]);
        blur(img, img, Size(3,3));
        img.convertTo(img, CV_32F, 1.0/255);

        Mat gray, lab, rgb, hsv, chn[3];
        cvtColor(img, img, COLOR_BGR2RGB);
        enhance(img, rgb);
        cvtColor(rgb, lab, COLOR_RGB2Lab);
        cvtColor(rgb, gray, COLOR_RGB2GRAY);
        cvtColor(rgb, hsv, COLOR_RGB2HSV);
        split(rgb,chn);


        Ptr<SuperpixelSLIC> slic = createSuperpixelSLIC(lab, algorithm, region_size, ruler);
        slic->iterate(num_iterations);
        if (min_element_size>0)
            slic->enforceLabelConnectivity(min_element_size);
        int n = slic->getNumberOfSuperpixels();
        Mat labels;
        slic->getLabels(labels);
        Scalar mg = mean(gray);
        std::vector<cluster> clusters;

        for (int j=0; j<n; j++) {
            Mat mask;
            compare(labels, j, mask, CMP_EQ);
            Moments mom = moments(mask);
            cluster c;
            c.label = 0;
            c.center = Point(mom.m10 / mom.m00, mom.m01 / mom.m00);
            Scalar m,d;
            meanStdDev(rgb, m,d, mask);
            c.features.push_back(m[0]);
            c.features.push_back(m[1]);
            c.features.push_back(m[2]);
            c.features.push_back(d[0]);
            c.features.push_back(d[1]);
            c.features.push_back(d[2]);
            meanStdDev(hsv, m,d, mask);
            c.features.push_back(m[0]);
            c.features.push_back(m[1]);
            c.features.push_back(d[1]);
            meanStdDev(gray, m,d, mask);
            c.features.push_back(m[0]);
            c.features.push_back(d[0]);
            double _m, _M;
            minMaxLoc(gray, &_m, &_M, 0,0, mask);
//            if (_m < 10)
//                continue;
            c.features.push_back(_m);
            c.features.push_back(_M);
            minMaxLoc(chn[0], &_m, &_M, 0,0, mask);
            c.features.push_back(_m);
            c.features.push_back(_M);
            minMaxLoc(chn[1], &_m, &_M, 0,0, mask);
            c.features.push_back(_m);
            c.features.push_back(_M);
            minMaxLoc(chn[2], &_m, &_M, 0,0, mask);
            c.features.push_back(_m);
            c.features.push_back(_M);
            clusters.push_back(c);
            //circle(img, center,2, Scalar(0,255,0),-1);
        }

        // normalize (whiten) features:
        std::vector<float> sums(20, 0), devs(20, 0);
        for (int j=0; j<n; j++) {
            double ctx = contextual(clusters, j, 120, mg[0]);
            clusters[j].features.push_back(ctx);
            for (int k=0; k<20; k++) {
                sums[k] += clusters[j].features.at<double>(k);
            }
        }
        for (int k=0; k<20; k++) {
            sums[k] /= clusters.size(); // mean
        }
        for (int j=0; j<n; j++) {
            for (int k=0; k<20; k++) {
                float diff = clusters[j].features.at<double>(k) - sums[k];
                devs[k] += diff*diff;
            }
        }
        for (int k=0; k<20; k++) {
            devs[k] /= clusters.size(); // stdev
            devs[k] = sqrt(devs[k]);
        }
        for (int j=0; j<n; j++) {
            for (int k=0; k<20; k++) {
                clusters[j].features.at<double>(k) -= sums[k];
                clusters[j].features.at<double>(k) /= devs[k];
            }
        }

        Point disc(-1,-1);
        // find the disc, and label the clusters
        for (int j=0; j<4; j++) {
            String fn(format("c:/data/diaret/gt/diaretdb1_image%03d_%02d_plain.xml",i+1,j+1));
            FileStorage xml(fn, 0);
            FileNode pnodes = xml.getFirstTopLevelNode(); //["markings"];
            for (FileNodeIterator it=pnodes.begin(); it!=pnodes.end(); ++it)
            {
                String type = (String)((*it)["type"]);
                String conf = (String)((*it)["conf"]);
                int x = (int)((*it)["x"]);
                int y = (int)((*it)["y"]);
                // lookup the cluster id from the slic labels obtained before
                int clusterid = labels.at<int>(y,x);
                if (type == "Disc")
                    disc = Point(x,y);
                if (type != "Hard_exudates")
                    continue;
                if (conf == "Low" || conf == "Min")
                    continue;
                clusters[clusterid].label = 1;
            }
        }

        fs << "clusters" << "[";
        for (int j=0; j<n; j++) {
            // cheating with the disc here ;)
            bool isdisc = norm(clusters[j].center - disc) < 110;
            fs << "{:";
            fs << "label" << (isdisc ? -1 : clusters[j].label);
            fs << "center" << clusters[j].center;
            fs << "features" << clusters[j].features;
            fs << "}";
        }
        fs << "]";
        fs << "}";

        cout << i << " " << n << " " << disc << endl;

        resize(rgb,rgb, Size(), .6f,.6f);
        cvtColor(rgb,rgb, COLOR_RGB2BGR);
        imshow("enhanced",rgb);

        // get the contours for displaying
        Mat mask;
        slic->getLabelContourMask(mask, true);
        img.setTo(Scalar(0, 0, 255), mask);

        resize(img, img, Size(), .6f,.6f);
        cvtColor(img,img, COLOR_RGB2BGR);
        imshow("Lab",img);
        if (waitKey(50)>0) break;
    }
    fs.release();
    return 0;
}

void prepare_data(Mat &trainData, Mat &trainLabels,
                  Mat &testData,  Mat &testLabels, std::vector<int> &testRanges) {
    std::ifstream in("c:/data/diaret/ddb1_v02_01_train_plain.txt");
    std::vector<String> trainNames;
    String line;
    while(getline(in,line)) {
        int imstop = line.find_first_of(' ');
        String img = line.substr(0,imstop);
        trainNames.push_back(img);
    }
    FileStorage fs("exudates.xml", 0);
    FileNode pnodes = fs.getFirstTopLevelNode();
    for (FileNodeIterator it=pnodes.begin(); it!=pnodes.end(); ++it)
    {
        String image = (String)((*it)["image"]);
        bool isTrain = false;
        for (size_t i=0; i<trainNames.size(); i++) {
            int istop1 = trainNames[i].find_last_of('/') + 1;
            int istop2 = image.find_last_of('\\') + 1;
            String s1 = trainNames[i].substr(istop1);
            String s2 = image.substr(istop2);
            if (s1 == s2) {
                isTrain = true;
                break;
            }
        }
        if (!isTrain)
            testRanges.push_back(0);

        FileNode clusters = (*it)["clusters"];
        for (FileNodeIterator cl=clusters.begin(); cl!=clusters.end(); ++cl)
        {
            int label = (int)((*cl)["label"]);
            if (label == -1) // ignore disc (this *should* be calculated at runtime)
                continue;
            Mat features;
            ((*cl)["features"]) >> features;
             if (isTrain) {
                 trainData.push_back(features.reshape(1,1));
                 trainLabels.push_back(label);
             } else {
                 testData.push_back(features.reshape(1,1));
                 testLabels.push_back(label);
                 testRanges.back() ++; // keep track, which image it belonged to
             }
         }
    }
    cout << trainData.size() << " " << testData.size() << endl;
}

void eval(const String &name, const Mat_<int> &confusion) {
    cout << confusion << endl;
    float TN = confusion(0,0);
    float TP = confusion(1,1);
    float FN = confusion(0,1);
    float FP = confusion(1,0);
    float acc = (TP+TN) / (TP+TN+FP+FN);
    float sen = TP / (TP+FN);
    float spc = TN / (TN+FP);
    cout << name << " " << " acc: " << acc << " sen: " << sen << " spc: " << spc << endl;
}

void test_svm() {
    std::vector<int> testRanges; // keep track, which image test features belonged to
    Mat trainData, trainLabels, testData, testLabels;
    prepare_data(trainData, trainLabels, testData, testLabels, testRanges);
    trainData.convertTo(trainData, CV_32F);
    testData.convertTo(testData, CV_32F);

    Ptr<ml::SVM> svm = ml::SVM::create();
    svm->setKernel(ml::SVM::INTER);
    svm->train(trainData, 0, trainLabels);
    int start=0;
    Mat_<int> confusion(2,2);
    for (size_t i=0; i<testRanges.size(); i++) {
        int end = testRanges[i]; // range per original image
        end += start;
        Mat td(testData, Range(start,end), Range::all());
        start = end;
        Mat res;
        svm->predict(td, res);
        for (int r=0; r<res.rows; r++) {
            int p = res.at<float>(r);
            int q = testLabels.at<int>(r);
            confusion(p,q) ++;
        }
    }
    eval("svm", confusion);
}

void test_lda() {
    std::vector<int> testRanges;
    Mat trainData, trainLabels, testData, testLabels;
    prepare_data(trainData, trainLabels, testData, testLabels, testRanges);
    LDA lda(trainData, trainLabels);
    Mat ev = lda.eigenvectors();
    int start=0;
    Mat_<int> confusion(2,2);
    for (size_t i=0; i<testRanges.size(); i++) {
        int end = testRanges[i];
        end += start;
        Mat td(testData, Range(start,end), Range::all());
        start = end;
        Mat proj = lda.project(td);
        for (int r=0; r<proj.rows; r++) {
            int p = proj.at<double>(r,0) > 0;
            int q = testLabels.at<int>(r);
            confusion(p,q) ++;
        }
    }
    eval("lda", confusion);
}

int main( int argc, char* argv[] )
{
    //make_dataset();
    //test_svm();
    test_lda();
    return 0;
}

to_fs.py

# process the groundtruth data, sowe can read it via opencv's FileStorage

import os,sys
import cv2, numpy as np
import xml.etree.ElementTree as ET # https://docs.python.org/3/library/xml.etree.elementtree.html


def convert(fn):
    print(fn)
    newname = fn.replace("groundtruth", "gt")
    fs = cv2.FileStorage(newname, 1)
    #help(cv2.FileNode)
    tree = ET.parse(fn)
    root = tree.getroot()
    fs.write("markings", "[")
    markings = root[1].findall("marking")
    print (len(markings))
    for m in markings:
        mt = m.find("markingtype").text
        cf = m.find("confidencelevel").text
        if m[0].tag =="polygonregion":
            coords = m[0].findall("coords2d")
        else:
            coords = m[0][0].findall("coords2d")
        print(m[0].tag, mt, len(coords), cf)
        for c in coords:
            cc = c.text.split(',')
            fs.write(None,"{:")
            fs.write("type",str(mt))
            fs.write("conf",str(cf))
            fs.write("x",int(cc[0]))
            fs.write("y",int(cc[1]))
            fs.write("}", None)
    fs.release()
    # cleanup
    f = open(newname, "rb")
    txt = f.read().decode(encoding='utf-8')
    txt = txt.replace("\"\"", "")
    txt = txt.replace("<type>", "<type>\"")
    txt = txt.replace("</type>", "\"</type>")
    txt = txt.replace("<conf>", "<conf>\"")
    txt = txt.replace("</conf>", "\"</conf>")
    f.close()
    f = open(newname,"wb")
    f.write(txt.encode(encoding='utf-8'))
    f.close()

def main():
    import glob
    xmls = glob.glob("c:/data/diaret/groundtruth/*.xml")
    for z in xmls:
        if z.find("_plain") > 0:
            convert(z)


main()

for the database groundtruth you only need to use images white and black no need to take the xml file you will find the image in : ...\diaretdb1_v_1_1\resources\images\ddb1_groundtruth\hardexudates
and could you explain please how the training work in your code i'm beginner and it's little hard to me to understand it ;)