Liczenie ziaren ryżu

81

Rozważ te 10 zdjęć różnych ilości niegotowanych ziaren białego ryżu.
TO TYLKO TRZEWIKI. Kliknij obraz, aby wyświetlić go w pełnym rozmiarze.

A: B: C: D: E:ZA b do re mi

F: G: H: I: J:fa sol H. ja jot

Ziarna: A: 3, B: 5, C: 12, D: 25, E: 50, F: 83, G: 120, H:150, I: 151, J: 200

Zauważ, że...

  • Ziarna mogą się stykać, ale nigdy się nie nakładają. Układ ziaren nigdy nie przekracza wysokości jednego ziarna.
  • Obrazy mają różne wymiary, ale skala ryżu we wszystkich jest spójna, ponieważ kamera i tło były nieruchome.
  • Ziarna nigdy nie wychodzą poza granice i nie dotykają granic obrazu.
  • Tło ma zawsze ten sam spójny odcień żółtawo-biały.
  • Małe i duże ziarna są liczone podobnie jak jedno ziarno.

Te 5 punktów jest gwarancją dla wszystkich obrazów tego rodzaju.

Wyzwanie

Napisz program, który pobiera takie obrazy i, tak dokładnie, jak to możliwe, liczy liczbę ziaren ryżu.

Twój program powinien pobrać nazwę pliku obrazu i wydrukować obliczoną liczbę ziaren. Twój program musi działać z co najmniej jednym z następujących formatów plików graficznych: JPEG, Bitmap, PNG, GIF, TIFF (obecnie wszystkie obrazy to JPEG).

Państwo może korzystać z bibliotek przetwarzania obrazu i wizyjne komputer.

Nie można zakodować na stałe wyników 10 przykładowych obrazów. Twój algorytm powinien mieć zastosowanie do wszystkich podobnych obrazów ziaren ryżu. Powinien być w stanie działać w mniej niż 5 minut na przyzwoitym nowoczesnym komputerze jeśli obszar obrazu jest mniejszy niż 2000 * 2000 pikseli i jest mniej niż 300 ziaren ryżu.

Punktacja

Dla każdego z 10 zdjęć weź wartość bezwzględną rzeczywistej liczby ziaren minus liczbę ziaren przewidzianych przez Twój program. Zsumuj te wartości bezwzględne, aby uzyskać wynik. Najniższy wynik wygrywa. Wynik 0 jest idealny.

W przypadku remisu wygrywa najwyższa głosowana odpowiedź. Mogę przetestować Twój program na dodatkowych obrazach, aby zweryfikować jego ważność i dokładność.

Hobby Calvina
źródło
1
Z pewnością ktoś musi spróbować scikit-learn!
Świetny konkurs! :) Btw - może nam powiedzieć coś o dacie zakończenia tego wyzwania?
cyriel,
1
@Lembik Down to 7 :)
Dr. belisarius
5
Pewnego dnia przyjedzie naukowiec zajmujący się ryżem i będzie szczęśliwy, że to pytanie istnieje.
Nit
2
@Nit Po prostu powiedz im ncbi.nlm.nih.gov/pmc/articles/PMC3510117 :)
Dr. belisarius

Odpowiedzi:

22

Mathematica, wynik: 7

i = {"http://i.stack.imgur.com/8T6W2.jpg",  "http://i.stack.imgur.com/pgWt1.jpg", 
     "http://i.stack.imgur.com/M0K5w.jpg",  "http://i.stack.imgur.com/eUFNo.jpg", 
     "http://i.stack.imgur.com/2TFdi.jpg",  "http://i.stack.imgur.com/wX48v.jpg", 
     "http://i.stack.imgur.com/eXCGt.jpg",  "http://i.stack.imgur.com/9na4J.jpg",
     "http://i.stack.imgur.com/UMP9V.jpg",  "http://i.stack.imgur.com/nP3Hr.jpg"};

im = Import /@ i;

Myślę, że nazwy funkcji są wystarczająco opisowe:

getSatHSVChannelAndBinarize[i_Image]             := Binarize@ColorSeparate[i, "HSB"][[2]]
removeSmallNoise[i_Image]                        := DeleteSmallComponents[i, 100]
fillSmallHoles[i_Image]                          := Closing[i, 1]
getMorphologicalComponentsAreas[i_Image]         := ComponentMeasurements[i, "Area"][[All, 2]]
roundAreaSizeToGrainCount[areaSize_, grainSize_] := Round[areaSize/grainSize]

Przetwarzanie wszystkich zdjęć jednocześnie:

counts = Plus @@@
  (roundAreaSizeToGrainCount[#, 2900] & /@
      (getMorphologicalComponentsAreas@
        fillSmallHoles@
         removeSmallNoise@
          getSatHSVChannelAndBinarize@#) & /@ im)

(* Output {3, 5, 12, 25, 49, 83, 118, 149, 152, 202} *)

Wynik to:

counts - {3, 5, 12, 25, 50, 83, 120, 150, 151, 200} // Abs // Total
(* 7 *)

Tutaj możesz zobaczyć czułość punktową względem użytej wielkości ziarna:

Grafika matematyczna

Dr Belizariusz
źródło
2
O wiele jaśniejsze, dzięki!
Calvin's Hobbies
Czy tę dokładną procedurę można skopiować w Pythonie, czy jest coś specjalnego w matematyce, czego nie mogą zrobić biblioteki Pythona?
@Lembik Nie mam pojęcia. Nie ma tutaj pytona. Przepraszam. (Jednak wątpię te same algorytmy EdgeDetect[], DeleteSmallComponents[]i Dilation[]są realizowane w innym miejscu)
dr Belizariusz
55

Python, wynik: 24 16

To rozwiązanie, podobnie jak Falko, polega na pomiarze obszaru „pierwszego planu” i podzieleniu go przez średni obszar ziarna.

W rzeczywistości program próbuje wykryć tło, a nie pierwszy plan. Biorąc pod uwagę fakt, że ziarna ryżu nigdy nie dotykają granicy obrazu, program rozpoczyna się od wypełnienia zalewającego białego w lewym górnym rogu. Algorytm zalewania maluje sąsiednie piksele, jeśli różnica między ich jasnością a bieżącym pikselem mieści się w określonym progu, dostosowując się w ten sposób do stopniowej zmiany koloru tła. Pod koniec tego etapu obraz może wyglądać mniej więcej tak:

Rycina 1

Jak widać, wykrywanie tła robi całkiem dobrą robotę, ale pomija obszary, które są „uwięzione” między ziarnami. Obsługujemy te obszary, szacując jasność tła na każdym pikselu i dopasowując wszystkie równe lub jaśniejsze piksele. Oszacowanie działa w ten sposób: podczas etapu wypełniania zalewamy obliczamy średnią jasność tła dla każdego wiersza i każdej kolumny. Szacowana jasność tła dla każdego piksela jest średnią jasnością wiersza i kolumny dla tego piksela. To powoduje coś takiego:

Rysunek 2

EDYCJA: W końcu obszar każdego ciągłego obszaru pierwszego planu (tj. Niebiałego) jest podzielony przez średnią, wstępnie obliczoną powierzchnię ziarna, dając nam oszacowanie liczby ziaren w tym regionie. Suma tych ilości jest wynikiem. Początkowo robiliśmy to samo dla całego obszaru pierwszego planu jako całości, ale to podejście jest dosłownie bardziej szczegółowe.


from sys import argv; from PIL import Image

# Init
I = Image.open(argv[1]); W, H = I.size; A = W * H
D = [sum(c) for c in I.getdata()]
Bh = [0] * H; Ch = [0] * H
Bv = [0] * W; Cv = [0] * W

# Flood-fill
Background = 3 * 255 + 1; S = [0]
while S:
    i = S.pop(); c = D[i]
    if c != Background:
        D[i] = Background
        Bh[i / W] += c; Ch[i / W] += 1
        Bv[i % W] += c; Cv[i % W] += 1
        S += [(i + o) % A for o in [1, -1, W, -W] if abs(D[(i + o) % A] - c) < 10]

# Eliminate "trapped" areas
for i in xrange(H): Bh[i] /= float(max(Ch[i], 1))
for i in xrange(W): Bv[i] /= float(max(Cv[i], 1))
for i in xrange(A):
    a = (Bh[i / W] + Bv[i % W]) / 2
    if D[i] >= a: D[i] = Background

# Estimate grain count
Foreground = -1; avg_grain_area = 3038.38; grain_count = 0
for i in xrange(A):
    if Foreground < D[i] < Background:
        S = [i]; area = 0
        while S:
            j = S.pop() % A
            if Foreground < D[j] < Background:
                D[j] = Foreground; area += 1
                S += [j - 1, j + 1, j - W, j + W]
        grain_count += int(round(area / avg_grain_area))

# Output
print grain_count

Wprowadza nazwę pliku wejściowego przez wiersz polecenia.

Wyniki

      Actual  Estimate  Abs. Error
A         3         3           0
B         5         5           0
C        12        12           0
D        25        25           0
E        50        48           2
F        83        83           0
G       120       116           4
H       150       145           5
I       151       156           5
J       200       200           0
                        ----------
                Total:         16

ZA b do re mi

fa sol H. ja jot

Łokieć
źródło
2
To naprawdę sprytne rozwiązanie, dobra robota!
Chris Cirefice,
1
skąd avg_grain_area = 3038.38;pochodzi?
njzk2
2
czy to się nie liczy hardcoding the result?
njzk2
5
@ njzk2 Nie. Biorąc pod uwagę regułę The images have different dimensions but the scale of the rice in all of them is consistent because the camera and background were stationary.Jest to tylko wartość reprezentująca tę regułę. Wynik zmienia się jednak zgodnie z danymi wejściowymi. Jeśli zmienisz regułę, wartość ta ulegnie zmianie, ale wynik będzie taki sam - na podstawie danych wejściowych.
Adam Davis
6
Nic mi nie jest ze średnim obszarem. Obszar ziarna jest (z grubsza) stały na wszystkich obrazach.
Calvin's Hobbies
28

Python + OpenCV: wynik 27

Skanowanie linii poziomej

Pomysł: zeskanuj obraz, jeden rząd na raz. Dla każdej linii policz liczbę napotkanych ziaren ryżu (sprawdzając, czy piksel zmienia kolor z czarnego na biały lub odwrotnie). Jeśli liczba ziaren dla linii wzrośnie (w porównaniu do poprzedniej linii), oznacza to, że napotkaliśmy nowe ziarno. Jeśli ta liczba spadnie, oznacza to, że przeszliśmy przez ziarno. W takim przypadku dodaj +1 do całkowitego wyniku.

wprowadź opis zdjęcia tutaj

Number in red = rice grains encountered for that line
Number in gray = total amount of grains encountered (what we are looking for)

Ze względu na sposób działania algorytmu ważne jest, aby mieć czysty, czarno-biały obraz. Dużo hałasu powoduje złe wyniki. Pierwsze główne tło jest czyszczone za pomocą wypełnienia zalewowego (rozwiązanie podobne do odpowiedzi Ell), a następnie stosuje się próg, aby uzyskać czarno-biały wynik.

wprowadź opis zdjęcia tutaj

Jest daleki od ideału, ale daje dobre wyniki w zakresie prostoty. Prawdopodobnie istnieje wiele sposobów na jego poprawę (poprzez zapewnienie lepszego obrazu czarno-białego, skanowanie w innych kierunkach (np .: pionowy, ukośny) przy średniej itp.)

import cv2
import numpy
import sys

filename = sys.argv[1]
I = cv2.imread(filename, 0)
h,w = I.shape[:2]
diff = (3,3,3)
mask = numpy.zeros((h+2,w+2),numpy.uint8)
cv2.floodFill(I,mask,(0,0), (255,255,255),diff,diff)
T,I = cv2.threshold(I,180,255,cv2.THRESH_BINARY)
I = cv2.medianBlur(I, 7)

totalrice = 0
oldlinecount = 0
for y in range(0, h):
    oldc = 0
    linecount = 0
    start = 0   
    for x in range(0, w):
        c = I[y,x] < 128;
        if c == 1 and oldc == 0:
            start = x
        if c == 0 and oldc == 1 and (x - start) > 10:
            linecount += 1
        oldc = c
    if oldlinecount != linecount:
        if linecount < oldlinecount:
            totalrice += oldlinecount - linecount
        oldlinecount = linecount
print totalrice

Błędy na obraz: 0, 0, 0, 3, 0, 12, 4, 0, 7, 1

tigrou
źródło
24

Python + OpenCV: Wynik 84

Oto pierwsza naiwna próba. Stosuje próg adaptacyjny z ręcznie dostrojonymi parametrami, zamyka niektóre otwory z późniejszą erozją i rozcieńczaniem oraz wyprowadza liczbę ziaren z obszaru pierwszego planu.

import cv2
import numpy as np

filename = raw_input()

I = cv2.imread(filename, 0)
I = cv2.medianBlur(I, 3)
bw = cv2.adaptiveThreshold(I, 255, cv2.ADAPTIVE_THRESH_GAUSSIAN_C, cv2.THRESH_BINARY, 101, 1)

kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (17, 17))
bw = cv2.dilate(cv2.erode(bw, kernel), kernel)

print np.round_(np.sum(bw == 0) / 3015.0)

Tutaj możesz zobaczyć pośrednie obrazy binarne (czarny jest na pierwszym planie):

wprowadź opis zdjęcia tutaj

Błędy na obraz to 0, 0, 2, 2, 4, 0, 27, 42, 0 i 7 ziaren.

Falko
źródło
20

C # + OpenCvSharp, wynik: 2

To moja druga próba. Zupełnie różni się od mojej pierwszej próby , która jest o wiele prostsza, więc zamieszczam to jako osobne rozwiązanie.

Podstawową ideą jest identyfikacja i oznakowanie każdego pojedynczego ziarna za pomocą iteracyjnego dopasowania elipsy. Następnie usuń piksele tego ziarna ze źródła i spróbuj znaleźć następne ziarno, aż każdy piksel zostanie oznaczony.

To nie jest najpiękniejsze rozwiązanie. To gigantyczny wieprz z 600 liniami kodu. Największy obraz potrzebuje 1,5 minuty. I naprawdę przepraszam za ten niechlujny kod.

Jest tak wiele parametrów i sposobów myślenia w tej rzeczy, że obawiam się, że mój program nie pasuje do 10 przykładowych obrazów. Końcowy wynik 2 jest prawie na pewno przypadkiem przeregulowania: mam dwa parametry average grain size in pixel, iminimum ratio of pixel / elipse_area , a na koniec po prostu wyczerpały wszystkie kombinacje tych dwóch parametrów, dopóki nie dostał najniższą ocenę. Nie jestem pewien, czy to wszystko koszerne z zasadami tego wyzwania.

average_grain_size_in_pixel = 2530
pixel / elipse_area >= 0.73

Ale nawet bez tych zbyt mocnych sprzęgieł wyniki są całkiem niezłe. Bez ustalonego rozmiaru ziarna lub stosunku pikseli, po prostu poprzez oszacowanie średniego rozmiaru ziarna na podstawie obrazów treningowych, wynik nadal wynosi 27.

I otrzymuję jako wynik nie tylko liczbę, ale rzeczywistą pozycję, orientację i kształt każdego ziarna. istnieje niewielka liczba źle oznakowanych ziaren, ale ogólnie większość etykiet dokładnie odpowiada rzeczywistym ziarnom:

A ZA B b C do D re Emi

F fa G sol H H. I ja Jjot

(kliknij na każde zdjęcie, aby wyświetlić wersję w pełnym rozmiarze)

Po tym etapie etykietowania mój program analizuje każde ziarno i ocenia na podstawie liczby pikseli i stosunku pikseli do pola elipsy, niezależnie od tego, czy jest to

  • pojedyncze ziarno (+1)
  • wiele ziaren błędnie oznakowanych jako jedno (+ X)
  • kropelka zbyt mała, aby być ziarnem (+0)

Oceny błędów dla każdego obrazu wynoszą A:0; B:0; C:0; D:0; E:2; F:0; G:0 ; H:0; I:0, J:0

Jednak rzeczywisty błąd jest prawdopodobnie nieco wyższy. Niektóre błędy w tym samym obrazie znoszą się nawzajem. W szczególności obraz H ma pewne źle oznakowane ziarna, podczas gdy na zdjęciu E etykiety są w większości prawidłowe

Pomysł jest trochę wymyślony:

  • Najpierw pierwszy plan jest oddzielany za pomocą progu otsu na kanale nasycenia (szczegóły w mojej poprzedniej odpowiedzi)

  • powtarzaj, aż nie pozostanie więcej pikseli:

    • wybierz największą kroplę
    • wybierz 10 losowych pikseli krawędzi na tej kropli jako pozycje początkowe dla ziarna

    • dla każdego punktu początkowego

      • załóż ziarno o wysokości i szerokości 10 pikseli w tej pozycji.

      • powtarzaj aż do konwergencji

        • idź od tego punktu promieniowo na zewnątrz, pod różnymi kątami, aż napotkasz piksel krawędzi (biały do ​​czarnego)

        • znalezione piksele powinny mieć nadzieję, że będą pikselami brzegowymi pojedynczego ziarna. Spróbuj oddzielić wartości odstające od wartości odstających, odrzucając piksele, które są bardziej odległe od zakładanej elipsy niż inne

        • wielokrotnie próbuj dopasować elipsę przez podzbiór wartości wewnętrznych, zachowaj najlepszą elipsę (RANSACK)

        • zaktualizuj położenie, orientację, szerokość i wysokość ziarna za pomocą znalezionej elipsy

        • jeśli pozycja ziarna nie zmienia się znacząco, zatrzymać

    • spośród 10 dopasowanych ziaren wybierz najlepsze ziarno według kształtu, liczby pikseli krawędzi. Odrzuć pozostałe

    • usuń wszystkie piksele tego ziarna z obrazu źródłowego, a następnie powtórz

    • na koniec przejrzyj listę znalezionych ziaren i policz każde ziarno jako 1 ziarno, 0 ziaren (za małe) lub 2 ziarna (za duże)

Jednym z moich głównych problemów było to, że nie chciałem implementować pełnej miary odległości w punkcie elipsy, ponieważ obliczenie, że samo w sobie jest skomplikowanym procesem iteracyjnym. Użyłem więc różnych obejść przy użyciu funkcji OpenCV Ellipse2Poly i FitEllipse, a wyniki nie są zbyt ładne.

Najwyraźniej przekroczyłem również limit rozmiaru dla codegolf.

Odpowiedź jest ograniczona do 30000 znaków, obecnie mam 34000. Więc będę musiał nieco skrócić poniższy kod.

Pełny kod można zobaczyć na stronie http://pastebin.com/RgM7hMxq

Przepraszam za to, nie wiedziałem, że istnieje limit rozmiaru.

class Program
{
    static void Main(string[] args)
    {

                // Due to size constraints, I removed the inital part of my program that does background separation. For the full source, check the link, or see my previous program.


                // list of recognized grains
                List<Grain> grains = new List<Grain>();

                Random rand = new Random(4); // determined by fair dice throw, guaranteed to be random

                // repeat until we have found all grains (to a maximum of 10000)
                for (int numIterations = 0; numIterations < 10000; numIterations++ )
                {
                    // erode the image of the remaining foreground pixels, only big blobs can be grains
                    foreground.Erode(erodedForeground,null,7);

                    // pick a number of starting points to fit grains
                    List<CvPoint> startPoints = new List<CvPoint>();
                    using (CvMemStorage storage = new CvMemStorage())
                    using (CvContourScanner scanner = new CvContourScanner(erodedForeground, storage, CvContour.SizeOf, ContourRetrieval.List, ContourChain.ApproxNone))
                    {
                        if (!scanner.Any()) break; // no grains left, finished!

                        // search for grains within the biggest blob first (this is arbitrary)
                        var biggestBlob = scanner.OrderByDescending(c => c.Count()).First();

                        // pick 10 random edge pixels
                        for (int i = 0; i < 10; i++)
                        {
                            startPoints.Add(biggestBlob.ElementAt(rand.Next(biggestBlob.Count())).Value);
                        }
                    }

                    // for each starting point, try to fit a grain there
                    ConcurrentBag<Grain> candidates = new ConcurrentBag<Grain>();
                    Parallel.ForEach(startPoints, point =>
                    {
                        Grain candidate = new Grain(point);
                        candidate.Fit(foreground);
                        candidates.Add(candidate);
                    });

                    Grain grain = candidates
                        .OrderByDescending(g=>g.Converged) // we don't want grains where the iterative fit did not finish
                        .ThenBy(g=>g.IsTooSmall) // we don't want tiny grains
                        .ThenByDescending(g => g.CircumferenceRatio) // we want grains that have many edge pixels close to the fitted elipse
                        .ThenBy(g => g.MeanSquaredError)
                        .First(); // we only want the best fit among the 10 candidates

                    // count the number of foreground pixels this grain has
                    grain.CountPixel(foreground);

                    // remove the grain from the foreground
                    grain.Draw(foreground,CvColor.Black);

                    // add the grain to the colection fo found grains
                    grains.Add(grain);
                    grain.Index = grains.Count;

                    // draw the grain for visualisation
                    grain.Draw(display, CvColor.Random());
                    grain.DrawContour(display, CvColor.Random());
                    grain.DrawEllipse(display, CvColor.Random());

                    //display.SaveImage("10-foundGrains.png");
                }

                // throw away really bad grains
                grains = grains.Where(g => g.PixelRatio >= 0.73).ToList();

                // estimate the average grain size, ignoring outliers
                double avgGrainSize =
                    grains.OrderBy(g => g.NumPixel).Skip(grains.Count/10).Take(grains.Count*9/10).Average(g => g.NumPixel);

                //ignore the estimated grain size, use a fixed size
                avgGrainSize = 2530;

                // count the number of grains, using the average grain size
                double numGrains = grains.Sum(g => Math.Round(g.NumPixel * 1.0 / avgGrainSize));

                // get some statistics
                double avgWidth = grains.Where(g => Math.Round(g.NumPixel * 1.0 / avgGrainSize) == 1).Average(g => g.Width);
                double avgHeight = grains.Where(g => Math.Round(g.NumPixel * 1.0 / avgGrainSize) == 1).Average(g => g.Height);
                double avgPixelRatio = grains.Where(g => Math.Round(g.NumPixel * 1.0 / avgGrainSize) == 1).Average(g => g.PixelRatio);

                int numUndersized = grains.Count(g => Math.Round(g.NumPixel * 1.0 / avgGrainSize) < 1);
                int numOversized = grains.Count(g => Math.Round(g.NumPixel * 1.0 / avgGrainSize) > 1);

                double avgWidthUndersized = grains.Where(g => Math.Round(g.NumPixel * 1.0 / avgGrainSize) < 1).Select(g=>g.Width).DefaultIfEmpty(0).Average();
                double avgHeightUndersized = grains.Where(g => Math.Round(g.NumPixel * 1.0 / avgGrainSize) < 1).Select(g => g.Height).DefaultIfEmpty(0).Average();
                double avgGrainSizeUndersized = grains.Where(g => Math.Round(g.NumPixel * 1.0 / avgGrainSize) < 1).Select(g => g.NumPixel).DefaultIfEmpty(0).Average();
                double avgPixelRatioUndersized = grains.Where(g => Math.Round(g.NumPixel * 1.0 / avgGrainSize) < 1).Select(g => g.PixelRatio).DefaultIfEmpty(0).Average();

                double avgWidthOversized = grains.Where(g => Math.Round(g.NumPixel * 1.0 / avgGrainSize) > 1).Select(g => g.Width).DefaultIfEmpty(0).Average();
                double avgHeightOversized = grains.Where(g => Math.Round(g.NumPixel * 1.0 / avgGrainSize) > 1).Select(g => g.Height).DefaultIfEmpty(0).Average();
                double avgGrainSizeOversized = grains.Where(g => Math.Round(g.NumPixel * 1.0 / avgGrainSize) > 1).Select(g => g.NumPixel).DefaultIfEmpty(0).Average();
                double avgPixelRatioOversized = grains.Where(g => Math.Round(g.NumPixel * 1.0 / avgGrainSize) > 1).Select(g => g.PixelRatio).DefaultIfEmpty(0).Average();


                Console.WriteLine("===============================");
                Console.WriteLine("Grains: {0}|{1:0.} of {2} (e{3}), size {4:0.}px, {5:0.}x{6:0.}  {7:0.000}  undersized:{8}  oversized:{9}   {10:0.0} minutes  {11:0.0} s per grain",grains.Count,numGrains,expectedGrains[fileNo],expectedGrains[fileNo]-numGrains,avgGrainSize,avgWidth,avgHeight, avgPixelRatio,numUndersized,numOversized,watch.Elapsed.TotalMinutes, watch.Elapsed.TotalSeconds/grains.Count);



                // draw the description for each grain
                foreach (Grain grain in grains)
                {
                    grain.DrawText(avgGrainSize, display, CvColor.Black);
                }

                display.SaveImage("10-foundGrains.png");
                display.SaveImage("X-" + file + "-foundgrains.png");
            }
        }
    }
}



public class Grain
{
    private const int MIN_WIDTH = 70;
    private const int MAX_WIDTH = 130;
    private const int MIN_HEIGHT = 20;
    private const int MAX_HEIGHT = 35;

    private static CvFont font01 = new CvFont(FontFace.HersheyPlain, 0.5, 1);
    private Random random = new Random(4); // determined by fair dice throw; guaranteed to be random


    /// <summary> center of grain </summary>
    public CvPoint2D32f Position { get; private set; }
    /// <summary> Width of grain (always bigger than height)</summary>
    public float Width { get; private set; }
    /// <summary> Height of grain (always smaller than width)</summary>
    public float Height { get; private set; }

    public float MinorRadius { get { return this.Height / 2; } }
    public float MajorRadius { get { return this.Width / 2; } }
    public double Angle { get; private set; }
    public double AngleRad { get { return this.Angle * Math.PI / 180; } }

    public int Index { get; set; }
    public bool Converged { get; private set; }
    public int NumIterations { get; private set; }
    public double CircumferenceRatio { get; private set; }
    public int NumPixel { get; private set; }
    public List<EllipsePoint> EdgePoints { get; private set; }
    public double MeanSquaredError { get; private set; }
    public double PixelRatio { get { return this.NumPixel / (Math.PI * this.MajorRadius * this.MinorRadius); } }
    public bool IsTooSmall { get { return this.Width < MIN_WIDTH || this.Height < MIN_HEIGHT; } }

    public Grain(CvPoint2D32f position)
    {
        this.Position = position;
        this.Angle = 0;
        this.Width = 10;
        this.Height = 10;
        this.MeanSquaredError = double.MaxValue;
    }

    /// <summary>  fit a single rice grain of elipsoid shape </summary>
    public void Fit(CvMat img)
    {
        // distance between the sampled points on the elipse circumference in degree
        int angularResolution = 1;

        // how many times did the fitted ellipse not change significantly?
        int numConverged = 0;

        // number of iterations for this fit
        int numIterations;

        // repeat until the fitted ellipse does not change anymore, or the maximum number of iterations is reached
        for (numIterations = 0; numIterations < 100 && !this.Converged; numIterations++)
        {
            // points on an ideal ellipse
            CvPoint[] points;
            Cv.Ellipse2Poly(this.Position, new CvSize2D32f(MajorRadius, MinorRadius), Convert.ToInt32(this.Angle), 0, 359, out points,
                            angularResolution);

            // points on the edge of foregroudn to background, that are close to the elipse
            CvPoint?[] edgePoints = new CvPoint?[points.Length];

            // remeber if the previous pixel in a given direction was foreground or background
            bool[] prevPixelWasForeground = new bool[points.Length];

            // when the first edge pixel is found, this value is updated
            double firstEdgePixelOffset = 200;

            // from the center of the elipse towards the outside:
            for (float offset = -this.MajorRadius + 1; offset < firstEdgePixelOffset + 20; offset++)
            {
                // draw an ellipse with the given offset
                Cv.Ellipse2Poly(this.Position, new CvSize2D32f(MajorRadius + offset, MinorRadius + (offset > 0 ? offset : MinorRadius / MajorRadius * offset)), Convert.ToInt32(this.Angle), 0,
                                359, out points, angularResolution);

                // for each angle
                Parallel.For(0, points.Length, i =>
                {
                    if (edgePoints[i].HasValue) return; // edge for this angle already found

                    // check if the current pixel is foreground
                    bool foreground = points[i].X < 0 || points[i].Y < 0 || points[i].X >= img.Cols || points[i].Y >= img.Rows
                                          ? false // pixel outside of image borders is always background
                                          : img.Get2D(points[i].Y, points[i].X).Val0 > 0;


                    if (prevPixelWasForeground[i] && !foreground)
                    {
                        // found edge pixel!
                        edgePoints[i] = points[i];

                        // if this is the first edge pixel we found, remember its offset. the other pixels cannot be too far away, so we can stop searching soon
                        if (offset < firstEdgePixelOffset && offset > 0) firstEdgePixelOffset = offset;
                    }

                    prevPixelWasForeground[i] = foreground;
                });
            }

            // estimate the distance of each found edge pixel from the ideal elipse
            // this is a hack, since the actual equations for estimating point-ellipse distnaces are complicated
            Cv.Ellipse2Poly(this.Position, new CvSize2D32f(MajorRadius, MinorRadius), Convert.ToInt32(this.Angle), 0, 360,
                            out points, angularResolution);
            var pointswithDistance =
                edgePoints.Select((p, i) => p.HasValue ? new EllipsePoint(p.Value, points[i], this.Position) : null)
                          .Where(p => p != null).ToList();

            if (pointswithDistance.Count == 0)
            {
                Console.WriteLine("no points found! should never happen! ");
                break;
            }

            // throw away all outliers that are too far outside the current ellipse
            double medianSignedDistance = pointswithDistance.OrderBy(p => p.SignedDistance).ElementAt(pointswithDistance.Count / 2).SignedDistance;
            var goodPoints = pointswithDistance.Where(p => p.SignedDistance < medianSignedDistance + 15).ToList();

            // do a sort of ransack fit with the inlier points to find a new better ellipse
            CvBox2D bestfit = ellipseRansack(goodPoints);

            // check if the fit has converged
            if (Math.Abs(this.Angle - bestfit.Angle) < 3 && // angle has not changed much (<3°)
                Math.Abs(this.Position.X - bestfit.Center.X) < 3 && // position has not changed much (<3 pixel)
                Math.Abs(this.Position.Y - bestfit.Center.Y) < 3)
            {
                numConverged++;
            }
            else
            {
                numConverged = 0;
            }

            if (numConverged > 2)
            {
                this.Converged = true;
            }

            //Console.WriteLine("Iteration {0}, delta {1:0.000} {2:0.000} {3:0.000}    {4:0.000}-{5:0.000} {6:0.000}-{7:0.000} {8:0.000}-{9:0.000}",
            //  numIterations, Math.Abs(this.Angle - bestfit.Angle), Math.Abs(this.Position.X - bestfit.Center.X), Math.Abs(this.Position.Y - bestfit.Center.Y), this.Angle, bestfit.Angle, this.Position.X, bestfit.Center.X, this.Position.Y, bestfit.Center.Y);

            double msr = goodPoints.Sum(p => p.Distance * p.Distance) / goodPoints.Count;

            // for drawing the polygon, filter the edge points more strongly
            if (goodPoints.Count(p => p.SignedDistance < 5) > goodPoints.Count / 2)
                goodPoints = goodPoints.Where(p => p.SignedDistance < 5).ToList();
            double cutoff = goodPoints.Select(p => p.Distance).OrderBy(d => d).ElementAt(goodPoints.Count * 9 / 10);
            goodPoints = goodPoints.Where(p => p.SignedDistance <= cutoff + 1).ToList();

            int numCertainEdgePoints = goodPoints.Count(p => p.SignedDistance > -2);
            this.CircumferenceRatio = numCertainEdgePoints * 1.0 / points.Count();

            this.Angle = bestfit.Angle;
            this.Position = bestfit.Center;
            this.Width = bestfit.Size.Width;
            this.Height = bestfit.Size.Height;
            this.EdgePoints = goodPoints;
            this.MeanSquaredError = msr;

        }
        this.NumIterations = numIterations;
        //Console.WriteLine("Grain found after {0,3} iterations, size={1,3:0.}x{2,3:0.}   pixel={3,5}    edgePoints={4,3}   msr={5,2:0.00000}", numIterations, this.Width,
        //                        this.Height, this.NumPixel, this.EdgePoints.Count, this.MeanSquaredError);
    }

    /// <summary> a sort of ransakc fit to find the best ellipse for the given points </summary>
    private CvBox2D ellipseRansack(List<EllipsePoint> points)
    {
        using (CvMemStorage storage = new CvMemStorage(0))
        {
            // calculate minimum bounding rectangle
            CvSeq<CvPoint> fullPointSeq = CvSeq<CvPoint>.FromArray(points.Select(p => p.Point), SeqType.EltypePoint, storage);
            var boundingRect = fullPointSeq.MinAreaRect2();

            // the initial candidate is the previously found ellipse
            CvBox2D bestEllipse = new CvBox2D(this.Position, new CvSize2D32f(this.Width, this.Height), (float)this.Angle);
            double bestError = calculateEllipseError(points, bestEllipse);

            Queue<EllipsePoint> permutation = new Queue<EllipsePoint>();
            if (points.Count >= 5) for (int i = -2; i < 20; i++)
                {
                    CvBox2D ellipse;
                    if (i == -2)
                    {
                        // first, try the ellipse described by the boundingg rect
                        ellipse = boundingRect;
                    }
                    else if (i == -1)
                    {
                        // then, try the best-fit ellipsethrough all points
                        ellipse = fullPointSeq.FitEllipse2();
                    }
                    else
                    {
                        // then, repeatedly fit an ellipse through a random sample of points

                        // pick some random points
                        if (permutation.Count < 5) permutation = new Queue<EllipsePoint>(permutation.Concat(points.OrderBy(p => random.Next())));
                        CvSeq<CvPoint> pointSeq = CvSeq<CvPoint>.FromArray(permutation.Take(10).Select(p => p.Point), SeqType.EltypePoint, storage);
                        for (int j = 0; j < pointSeq.Count(); j++) permutation.Dequeue();

                        // fit an ellipse through these points
                        ellipse = pointSeq.FitEllipse2();
                    }

                    // assure that the width is greater than the height
                    ellipse = NormalizeEllipse(ellipse);

                    // if the ellipse is too big for agrain, shrink it
                    ellipse = rightSize(ellipse, points.Where(p => isOnEllipse(p.Point, ellipse, 10, 10)).ToList());

                    // sometimes the ellipse given by FitEllipse2 is totally off
                    if (boundingRect.Center.DistanceTo(ellipse.Center) > Math.Max(boundingRect.Size.Width, boundingRect.Size.Height) * 2)
                    {
                        // ignore this bad fit
                        continue;
                    }

                    // estimate the error
                    double error = calculateEllipseError(points, ellipse);

                    if (error < bestError)
                    {
                        // found a better ellipse!
                        bestError = error;
                        bestEllipse = ellipse;
                    }
                }

            return bestEllipse;
        }
    }

    /// <summary> The proper thing to do would be to use the actual distance of each point to the elipse.
    /// However that formula is complicated, so ...  </summary>
    private double calculateEllipseError(List<EllipsePoint> points, CvBox2D ellipse)
    {
        const double toleranceInner = 5;
        const double toleranceOuter = 10;
        int numWrongPoints = points.Count(p => !isOnEllipse(p.Point, ellipse, toleranceInner, toleranceOuter));
        double ratioWrongPoints = numWrongPoints * 1.0 / points.Count;

        int numTotallyWrongPoints = points.Count(p => !isOnEllipse(p.Point, ellipse, 10, 20));
        double ratioTotallyWrongPoints = numTotallyWrongPoints * 1.0 / points.Count;

        // this pseudo-distance is biased towards deviations on the major axis
        double pseudoDistance = Math.Sqrt(points.Sum(p => Math.Abs(1 - ellipseMetric(p.Point, ellipse))) / points.Count);

        // primarily take the number of points far from the elipse border as an error metric.
        // use pseudo-distance to break ties between elipses with the same number of wrong points
        return ratioWrongPoints * 1000  + ratioTotallyWrongPoints+ pseudoDistance / 1000;
    }


    /// <summary> shrink an ellipse if it is larger than the maximum grain dimensions </summary>
    private static CvBox2D rightSize(CvBox2D ellipse, List<EllipsePoint> points)
    {
        if (ellipse.Size.Width < MAX_WIDTH && ellipse.Size.Height < MAX_HEIGHT) return ellipse;

        // elipse is bigger than the maximum grain size
        // resize it so it fits, while keeping one edge of the bounding rectangle constant

        double desiredWidth = Math.Max(10, Math.Min(MAX_WIDTH, ellipse.Size.Width));
        double desiredHeight = Math.Max(10, Math.Min(MAX_HEIGHT, ellipse.Size.Height));

        CvPoint2D32f average = points.Average();

        // get the corners of the surrounding bounding box
        var corners = ellipse.BoxPoints().ToList();

        // find the corner that is closest to the center of mass of the points
        int i0 = ellipse.BoxPoints().Select((point, index) => new { point, index }).OrderBy(p => p.point.DistanceTo(average)).First().index;
        CvPoint p0 = corners[i0];

        // find the two corners that are neighbouring this one
        CvPoint p1 = corners[(i0 + 1) % 4];
        CvPoint p2 = corners[(i0 + 3) % 4];

        // p1 is the next corner along the major axis (widht), p2 is the next corner along the minor axis (height)
        if (p0.DistanceTo(p1) < p0.DistanceTo(p2))
        {
            CvPoint swap = p1;
            p1 = p2;
            p2 = swap;
        }

        // calculate the three other corners with the desired widht and height

        CvPoint2D32f edge1 = (p1 - p0);
        CvPoint2D32f edge2 = p2 - p0;
        double edge1Length = Math.Max(0.0001, p0.DistanceTo(p1));
        double edge2Length = Math.Max(0.0001, p0.DistanceTo(p2));

        CvPoint2D32f newCenter = (CvPoint2D32f)p0 + edge1 * (desiredWidth / edge1Length) + edge2 * (desiredHeight / edge2Length);

        CvBox2D smallEllipse = new CvBox2D(newCenter, new CvSize2D32f((float)desiredWidth, (float)desiredHeight), ellipse.Angle);

        return smallEllipse;
    }

    /// <summary> assure that the width of the elipse is the major axis, and the height is the minor axis.
    /// Swap widht/height and rotate by 90° otherwise  </summary>
    private static CvBox2D NormalizeEllipse(CvBox2D ellipse)
    {
        if (ellipse.Size.Width < ellipse.Size.Height)
        {
            ellipse = new CvBox2D(ellipse.Center, new CvSize2D32f(ellipse.Size.Height, ellipse.Size.Width), (ellipse.Angle + 90 + 360) % 360);
        }
        return ellipse;
    }

    /// <summary> greater than 1 for points outside ellipse, smaller than 1 for points inside ellipse </summary>
    private static double ellipseMetric(CvPoint p, CvBox2D ellipse)
    {
        double theta = ellipse.Angle * Math.PI / 180;
        double u = Math.Cos(theta) * (p.X - ellipse.Center.X) + Math.Sin(theta) * (p.Y - ellipse.Center.Y);
        double v = -Math.Sin(theta) * (p.X - ellipse.Center.X) + Math.Cos(theta) * (p.Y - ellipse.Center.Y);

        return u * u / (ellipse.Size.Width * ellipse.Size.Width / 4) + v * v / (ellipse.Size.Height * ellipse.Size.Height / 4);
    }

    /// <summary> Is the point on the ellipseBorder, within a certain tolerance </summary>
    private static bool isOnEllipse(CvPoint p, CvBox2D ellipse, double toleranceInner, double toleranceOuter)
    {
        double theta = ellipse.Angle * Math.PI / 180;
        double u = Math.Cos(theta) * (p.X - ellipse.Center.X) + Math.Sin(theta) * (p.Y - ellipse.Center.Y);
        double v = -Math.Sin(theta) * (p.X - ellipse.Center.X) + Math.Cos(theta) * (p.Y - ellipse.Center.Y);

        double innerEllipseMajor = (ellipse.Size.Width - toleranceInner) / 2;
        double innerEllipseMinor = (ellipse.Size.Height - toleranceInner) / 2;
        double outerEllipseMajor = (ellipse.Size.Width + toleranceOuter) / 2;
        double outerEllipseMinor = (ellipse.Size.Height + toleranceOuter) / 2;

        double inside = u * u / (innerEllipseMajor * innerEllipseMajor) + v * v / (innerEllipseMinor * innerEllipseMinor);
        double outside = u * u / (outerEllipseMajor * outerEllipseMajor) + v * v / (outerEllipseMinor * outerEllipseMinor);
        return inside >= 1 && outside <= 1;
    }


    /// <summary> count the number of foreground pixels for this grain </summary>
    public int CountPixel(CvMat img)
    {
        // todo: this is an incredibly inefficient way to count, allocating a new image with the size of the input each time
        using (CvMat mask = new CvMat(img.Rows, img.Cols, MatrixType.U8C1))
        {
            mask.SetZero();
            mask.FillPoly(new CvPoint[][] { this.EdgePoints.Select(p => p.Point).ToArray() }, CvColor.White);
            mask.And(img, mask);
            this.NumPixel = mask.CountNonZero();
        }
        return this.NumPixel;
    }

    /// <summary> draw the recognized shape of the grain </summary>
    public void Draw(CvMat img, CvColor color)
    {
        img.FillPoly(new CvPoint[][] { this.EdgePoints.Select(p => p.Point).ToArray() }, color);
    }

    /// <summary> draw the contours of the grain </summary>
    public void DrawContour(CvMat img, CvColor color)
    {
        img.DrawPolyLine(new CvPoint[][] { this.EdgePoints.Select(p => p.Point).ToArray() }, true, color);
    }

    /// <summary> draw the best-fit ellipse of the grain </summary>
    public void DrawEllipse(CvMat img, CvColor color)
    {
        img.DrawEllipse(this.Position, new CvSize2D32f(this.MajorRadius, this.MinorRadius), this.Angle, 0, 360, color, 1);
    }

    /// <summary> print the grain index and the number of pixels divided by the average grain size</summary>
    public void DrawText(double averageGrainSize, CvMat img, CvColor color)
    {
        img.PutText(String.Format("{0}|{1:0.0}", this.Index, this.NumPixel / averageGrainSize), this.Position + new CvPoint2D32f(-5, 10), font01, color);
    }

}

Jestem trochę zawstydzony tym rozwiązaniem, ponieważ a) nie jestem pewien, czy jest to zgodne z duchem tego wyzwania, i b) jest zbyt duże, aby odpowiedzieć na codegolf i brakuje mu elegancji innych rozwiązań.

Z drugiej strony jestem całkiem zadowolony z postępu, jaki osiągnąłem w etykietowaniu ziaren, a nie tylko ich liczeniu, więc jest.

HugoRune
źródło
Wiesz, że możesz skrócić ten kod o wielkości, używając mniejszych nazw i innych technik gry w golfa;)
Optimizer
Prawdopodobnie, ale nie chciałem dalej zaciemniać tego rozwiązania. Jest to zbyt zaciemnione jak na mój gust :)
HugoRune
+1 za wysiłek i ponieważ jesteś jedynym, który znajduje sposób, aby wyświetlić indywidualnie każde ziarno. Niestety kod jest nieco rozdęty i polega na wielu na stałych zapisanych na stałe. Byłbym ciekawy, jak algorytm skanowania linii, który napisałem, działa na tym (na poszczególnych kolorowych ziarnach).
tigrou
Naprawdę uważam, że jest to właściwe podejście do tego rodzaju problemu (+1 dla ciebie), ale zastanawiam się, dlaczego „wybierasz 10 losowych pikseli krawędzi”, pomyślałbym, że uzyskałbyś lepszą wydajność, gdybyś wybrał punkty krawędzi z najmniejszą liczbą pobliskich punktów krawędzi (tj. części, które wystają), pomyślałem (teoretycznie) to wyeliminowałoby najpierw „najłatwiejsze” ziarna, czy zastanawiałeś się nad tym?
David Rogers
Myślałem o tym, ale jeszcze tego nie próbowałem. „10 losowych pozycji początkowych” było późnym dodatkiem, który łatwo było dodać i zrównoleglić. Wcześniej „jedna losowa pozycja początkowa” była znacznie lepsza niż „zawsze lewy górny róg”. Niebezpieczeństwo wyboru pozycji początkowych za pomocą tej samej strategii za każdym razem polega na tym, że kiedy usunę najlepsze dopasowanie, pozostałe 9 prawdopodobnie zostanie ponownie wybrane następnym razem, a wraz z upływem czasu najgorsze z tych pozycji początkowych pozostaną w tyle i zostaną ponownie wybrane i jeszcze raz. Częścią, która wystaje, mogą być tylko resztki niecałkowicie usuniętego poprzedniego ziarna.
HugoRune,
17

C ++, OpenCV, wynik: 9

Podstawowa idea mojej metody jest dość prosta - spróbuj usunąć pojedyncze ziarna (i „podwójne ziarna” - 2 (ale nie więcej!) Ziarna, blisko siebie) z obrazu, a następnie policz resztę za pomocą metody opartej na powierzchni (np. Falko, Ell i Belizariusz). Korzystanie z tego podejścia jest nieco lepsze niż standardowa „metoda obszarowa”, ponieważ łatwiej jest znaleźć dobrą średnią wartośćPixelsPerObject.

(1. krok) Przede wszystkim musimy zastosować binaryzację Otsu na kanale S obrazu w HSV. Kolejnym krokiem jest użycie operatora dylatacji w celu poprawy jakości wyodrębnionego pierwszego planu. Niż musimy znaleźć kontury. Oczywiście niektóre kontury nie są ziarnami ryżu - musimy usunąć zbyt małe kontury (o powierzchni mniejszej niż średniaPixelsPerObject / 4. W mojej sytuacji średniaPixelsPerObject to 2855). Teraz w końcu możemy zacząć zliczać ziarna :) (drugi krok) Znalezienie pojedynczych i podwójnych ziaren jest dość proste - wystarczy spojrzeć na liście konturów, aby uzyskać kontury o obszarze w określonych zakresach - jeśli obszar konturu jest w zasięgu, usuń go z listy i dodaj 1 (lub 2, jeśli było to „podwójne” ziarno) do licznika ziaren. (3. krok) Ostatnim krokiem jest oczywiście podzielenie obszaru pozostałych konturów przez wartość średniąPixelsPerObject i dodanie wyniku do licznika ziaren.

Obrazy (dla obrazu F.jpg) powinny lepiej przedstawiać ten pomysł niż słowa:
1. krok (bez małych konturów (szumów)): 1. stopień (bez małych konturów (hałas))
2. krok - tylko proste kontury: Drugi krok - tylko proste kontury
3. krok - pozostałe kontury: Trzeci krok - pozostałe kontury

Oto kod, jest raczej brzydki, ale powinien działać bez problemu. Oczywiście wymagany jest OpenCV.

#include "stdafx.h"

#include <cv.hpp>
#include <cxcore.h>
#include <highgui.h>
#include <vector>

using namespace cv;
using namespace std;

//A: 3, B: 5, C: 12, D: 25, E: 50, F: 83, G: 120, H:150, I: 151, J: 200
const int goodResults[] = {3, 5, 12, 25, 50, 83, 120, 150, 151, 200};
const float averagePixelsPerObject = 2855.0;

const int singleObjectPixelsCountMin = 2320;
const int singleObjectPixelsCountMax = 4060;

const int doubleObjectPixelsCountMin = 5000;
const int doubleObjectPixelsCountMax = 8000;

float round(float x)
{
    return x >= 0.0f ? floorf(x + 0.5f) : ceilf(x - 0.5f);
}

Mat processImage(Mat m, int imageIndex, int &error)
{
    int objectsCount = 0;
    Mat output, thresholded;
    cvtColor(m, output, CV_BGR2HSV);
    vector<Mat> channels;
    split(output, channels);
    threshold(channels[1], thresholded, 0, 255, CV_THRESH_OTSU | CV_THRESH_BINARY);
    dilate(thresholded, output, Mat()); //dilate to imporove quality of binary image
    imshow("thresholded", thresholded);
    int nonZero = countNonZero(output); //not realy important - just for tests
    if (imageIndex != -1)
        cout << "non zero: " << nonZero << ", average pixels per object: " << nonZero/goodResults[imageIndex] << endl;
    else
        cout << "non zero: " << nonZero << endl;

    vector<vector<Point>> contours, contoursOnlyBig, contoursWithoutSimpleObjects, contoursSimple;
    findContours(output, contours, CV_RETR_EXTERNAL, CV_CHAIN_APPROX_NONE); //find only external contours
    for (int i=0; i<contours.size(); i++)
        if (contourArea(contours[i]) > averagePixelsPerObject/4.0)
            contoursOnlyBig.push_back(contours[i]); //add only contours with area > averagePixelsPerObject/4 ---> skip small contours (noise)

    Mat bigContoursOnly = Mat::zeros(output.size(), output.type());
    Mat allContours = bigContoursOnly.clone();
    drawContours(allContours, contours, -1, CV_RGB(255, 255, 255), -1);
    drawContours(bigContoursOnly, contoursOnlyBig, -1, CV_RGB(255, 255, 255), -1);
    //imshow("all contours", allContours);
    output = bigContoursOnly;

    nonZero = countNonZero(output); //not realy important - just for tests
    if (imageIndex != -1)
        cout << "non zero: " << nonZero << ", average pixels per object: " << nonZero/goodResults[imageIndex] << " objects: "  << goodResults[imageIndex] << endl;
    else
        cout << "non zero: " << nonZero << endl;

    for (int i=0; i<contoursOnlyBig.size(); i++)
    {
        double area = contourArea(contoursOnlyBig[i]);
        if (area >= singleObjectPixelsCountMin && area <= singleObjectPixelsCountMax) //is this contours a single grain ?
        {
            contoursSimple.push_back(contoursOnlyBig[i]);
            objectsCount++;
        }
        else
        {
            if (area >= doubleObjectPixelsCountMin && area <= doubleObjectPixelsCountMax) //is this contours a double grain ?
            {
                contoursSimple.push_back(contoursOnlyBig[i]);
                objectsCount+=2;
            }
            else
                contoursWithoutSimpleObjects.push_back(contoursOnlyBig[i]); //group of grainss
        }
    }

    cout << "founded single objects: " << objectsCount << endl;
    Mat thresholdedImageMask = Mat::zeros(output.size(), output.type()), simpleContoursMat = Mat::zeros(output.size(), output.type());
    drawContours(simpleContoursMat, contoursSimple, -1, CV_RGB(255, 255, 255), -1);
    if (contoursWithoutSimpleObjects.size())
        drawContours(thresholdedImageMask, contoursWithoutSimpleObjects, -1, CV_RGB(255, 255, 255), -1); //draw only contours of groups of grains
    imshow("simpleContoursMat", simpleContoursMat);
    imshow("thresholded image mask", thresholdedImageMask);
    Mat finalResult;
    thresholded.copyTo(finalResult, thresholdedImageMask); //copy using mask - only pixels whc=ich belongs to groups of grains will be copied
    //imshow("finalResult", finalResult);
    nonZero = countNonZero(finalResult); // count number of pixels in all gropus of grains (of course without single or double grains)
    int goodObjectsLeft = goodResults[imageIndex]-objectsCount;
    if (imageIndex != -1)
        cout << "non zero: " << nonZero << ", average pixels per object: " << (goodObjectsLeft ? (nonZero/goodObjectsLeft) : 0) << " objects left: " << goodObjectsLeft <<  endl;
    else
        cout << "non zero: " << nonZero << endl;
    objectsCount += round((float)nonZero/(float)averagePixelsPerObject);

    if (imageIndex != -1)
    {
        error = objectsCount-goodResults[imageIndex];
        cout << "final objects count: " << objectsCount << ", should be: " << goodResults[imageIndex] << ", error is: " << error <<  endl;
    }
    else
        cout << "final objects count: " << objectsCount << endl; 
    return output;
}

int main(int argc, char* argv[])
{
    string fileName = "A";
    int totalError = 0, error;
    bool fastProcessing = true;
    vector<int> errors;

    if (argc > 1)
    {
        Mat m = imread(argv[1]);
        imshow("image", m);
        processImage(m, -1, error);
        waitKey(-1);
        return 0;
    }

    while(true)
    {
        Mat m = imread("images\\" + fileName + ".jpg");
        cout << "Processing image: " << fileName << endl;
        imshow("image", m);
        processImage(m, fileName[0] - 'A', error);
        totalError += abs(error);
        errors.push_back(error);
        if (!fastProcessing && waitKey(-1) == 'q')
            break;
        fileName[0] += 1;
        if (fileName[0] > 'J')
        {
            if (fastProcessing)
                break;
            else
                fileName[0] = 'A';
        }
    }
    cout << "Total error: " << totalError << endl;
    cout << "Errors: " << (Mat)errors << endl;
    cout << "averagePixelsPerObject:" << averagePixelsPerObject << endl;

    return 0;
}

Jeśli chcesz zobaczyć wyniki wszystkich kroków, odkomentuj wszystkie wywołania funkcji imshow (.., ..) i ustaw zmienną fastProcessing na wartość false. Obrazy (A.jpg, B.jpg, ...) powinny znajdować się w obrazach katalogowych. Alternatywnie możesz podać nazwę jednego obrazu jako parametr z wiersza poleceń.

Oczywiście, jeśli coś jest niejasne, mogę to wyjaśnić i / lub dostarczyć trochę zdjęć / informacji.

cyriel
źródło
12

C # + OpenCvSharp, wynik: 71

To jest bardzo irytujące, próbowałem znaleźć rozwiązanie, które faktycznie identyfikuje każde ziarno za pomocą zlewu , ale po prostu. żargon. otrzymać. to. do. praca.

Postanowiłem znaleźć rozwiązanie, które przynajmniej rozdzieli niektóre pojedyncze ziarna, a następnie wykorzystuje je do oszacowania średniej wielkości ziarna. Jednak do tej pory nie jestem w stanie pokonać rozwiązań o twardej ziarnistości.

Główną zaletą tego rozwiązania jest to, że nie zakłada ono stałego rozmiaru piksela dla ziaren i powinno działać nawet po przesunięciu aparatu lub zmianie rodzaju ryżu.

A.jpg; liczba ziaren: 3; oczekiwany 3; błąd 0; liczba pikseli na ziarno: 2525,0;
B.jpg; liczba ziaren: 7; oczekiwany 5; błąd 2; liczba pikseli na ziarno: 1920,0;
C.jpg; liczba ziaren: 6; spodziewany 12; błąd 6; liczba pikseli na ziarno: 4242,5;
D.jpg; liczba ziaren: 23; spodziewany 25; błąd 2; liczba pikseli na ziarno: 2415,5;
E.jpg; liczba ziaren: 47; spodziewane 50; błąd 3; liczba pikseli na ziarno: 2729,9;
F.jpg; liczba ziaren: 65; oczekiwany 83; błąd 18; liczba pikseli na ziarno: 2860,5;
G.jpg; liczba ziaren: 120; spodziewane 120; błąd 0; liczba pikseli na ziarno: 2552,3;
H.jpg; liczba ziaren: 159; spodziewany 150; błąd 9; liczba pikseli na ziarno: 2624,7;
I.jpg; liczba ziaren: 141; oczekiwany 151; błąd 10; liczba pikseli na ziarno: 2697,4;
J.jpg; liczba ziaren: 179; spodziewane 200; błąd 21; liczba pikseli na ziarno: 2847,1;
całkowity błąd: 71

Moje rozwiązanie działa w następujący sposób:

Rozdziel pierwszy plan, przekształcając obraz do HSV i stosując próg Otsu w kanale nasycenia. Jest to bardzo proste, działa bardzo dobrze i poleciłbym to wszystkim, którzy chcą spróbować tego wyzwania:

saturation channel                -->         Otsu thresholding

wprowadź opis zdjęcia tutaj -> wprowadź opis zdjęcia tutaj

Spowoduje to czyste usunięcie tła.

Następnie dodatkowo usunąłem cienie ziarna z pierwszego planu, stosując stały próg do kanału wartości. (Nie jestem pewien, czy to naprawdę pomaga, ale wystarczyło dodać)

wprowadź opis zdjęcia tutaj

Następnie stosuję transformację odległości na obrazie na pierwszym planie.

wprowadź opis zdjęcia tutaj

i znajdź wszystkie lokalne maksima w tej transformacji odległości.

To tutaj rozpada się mój pomysł. aby uniknąć mnożenia lokalnych maksimów w obrębie tego samego ziarna, muszę dużo filtrować. Obecnie utrzymuję tylko najsilniejsze maksimum w promieniu 45 pikseli, co oznacza, że ​​nie każde ziarno ma lokalne maksimum. I tak naprawdę nie mam uzasadnienia dla promienia 45 pikseli, to była tylko wartość, która działała.

wprowadź opis zdjęcia tutaj

(jak widać, nie jest to wystarczająco dużo nasion, aby uwzględnić każde ziarno)

Następnie używam tych maksimów jako nasion dla algorytmu zlewu:

wprowadź opis zdjęcia tutaj

Wyniki są meh . Miałem nadzieję na pojedyncze ziarna, ale kępy są wciąż zbyt duże.

Teraz identyfikuję najmniejsze plamy, liczę ich średni rozmiar w pikselach, a następnie szacuję z nich liczbę ziaren. To nie to, co planuje zrobić na początku, ale to był jedyny sposób, by uratować to.

using System ; 
using System . Kolekcje . Ogólne ; 
using System . Linq ; 
using System . Tekst ; 
using OpenCvSharp ;

przestrzeń nazw GrainTest2 { program klasy { static void Main ( string [] args ) { string [] files = new [] { "sourceA.jpg" , "sourceB.jpg" , "sourceC.jpg" , "sourceD.jpg" , " sourceE.jpg ” , „ sourceF.jpg ” , „ sourceG.jpg ” , „ sourceH.jpg ” , „ sourceI.jpg ” , „ sourceJ.jpg ” , };int [] expectGrains

     
    
          
        
             
                               
                                     
                                     
                                      
                               
            = nowy [] { 3 , 5 , 12 , 25 , 50 , 83 , 120 , 150 , 151 , 200 ,};          

            int totalError = 0 ; int totalPixels = 0 ; 
             

            for ( int fileNo = 0 ; fileNo markers = new List (); { // ustaw każde lokalne maksimum jako numer początkowy 25, 35, 45, ... // (liczby rzeczywiste nie mają znaczenia, wybrane dla lepszej widoczności w png ) int marker No    
                    using ( CvMemStorage storage = new CvMemStorage ()) 
                    using ( CvContourScanner scanner = new CvContourScanner ( localMaxima , storage , CvContour . SizeOf , ContourRetrieval . External , ContourChain . ApproxNone ))      
                    
                        
                        
                         = 20 ; foreach ( CvSeq c w skanerze ) { 
                            markerNo + = 5 ; 
                            markery . Dodaj ( markerNo ); 
                            waterShedMarkers . DrawContours ( c , nowy CvScalar ( markerNo ), nowy 
                         
                             CvScalar ( markerNo ), 0 , - 1 ); } } 
                    waterShedMarkers . SaveImage ( „08-watershed-seeds.png” );  
                        
                    


                    źródło . Watershed ( waterShedMarkers ); 
                    waterShedMarkers . SaveImage ( „09-watershed-result.png” );


                    List pixelPerBlob = new List ();  

                    // Okropny hack, ponieważ nie mogłem zmusić Cv2.ConnectedComponents do pracy z tym opakowaniem openCv // Więc podjąłem obejście, aby policzyć liczbę pikseli na 
                    kroplę waterShedMarkers . ConvertScale ( waterShedThreshold ); foreach ( Int markerNo w markerów ) { 
                        stosując ( { 
                            waterShedMarkers . CMP
                    
                     
                    CvMat tmp = new CvMat(waterShedMarkers.Rows, waterShedThreshold.Cols, MatrixType.U8C1))
                        (markerNo, tmp, ArrComparison.EQ);
                            pixelsPerBlob.Add(tmp.CountNonZero());

                        }
                    }

                    // oszacuj wielkość pojedynczego ziarna // krok 1: załóż, że 10% najmniejsza kropelka to całe ziarno; podwójny singleGrain
                    
                    = pixelsPerBlob.OrderBy(p => p).ElementAt(pixelsPerBlob.Count/15);

                    // step2: take all blobs that are not much bigger than the currently estimated singel grain size
                    //        average their size
                    //        repeat until convergence (too lazy to check for convergence)
                    for (int i = 0; i  p  Math.Round(p/singleGrain)).Sum());

                    Console.WriteLine("input: {0}; number of grains: {1,4:0.}; expected {2,4}; error {3,4}; pixels per grain: {4:0.0}; better: {5:0.}", file, numGrains, expectedGrains[fileNo], Math.Abs(numGrains - expectedGrains[fileNo]), singleGrain, pixelsPerBlob.Sum() / 1434.9);

                    totalError += Math.Abs(numGrains - expectedGrains[fileNo]);
                    totalPixels += pixelsPerBlob.Sum();

                    // this is a terrible hack to visualise the estimated number of grains per blob.
                    // i'm too tired to clean it up
                    #region please ignore
                    using (CvMemStorage storage = new CvMemStorage())
                    using (CvMat tmp = waterShedThreshold.Clone())
                    using (CvMat tmpvisu = new CvMat(source.Rows, source.Cols, MatrixType.S8C3))
                    {
                        foreach (int markerNo in markers)
                        {
                            tmp.SetZero();
                            waterShedMarkers.CmpS(markerNo, tmp, ArrComparison.EQ);
                            double curGrains = tmp.CountNonZero() * 1.0 / singleGrain;
                            using (
                                CvContourScanner scanner = new CvContourScanner(tmp, storage, CvContour.SizeOf, ContourRetrieval . Zewnętrzny , ContourChain .  
                                                                                ApproxNone))
                            {
                                tmpvisu.Set(CvColor.Random PutText ( „” + Math . Round ( curGrains , 1 ), c . First (). Wartość , nowy CvFont ( FontFace . HersheyPlain(), tmp);
                                foreach (CvSeq c in scanner)
                                {
                                    //tmpvisu.DrawContours(c, CvColor.Random(), CvColor.DarkGreen, 0, -1);
                                    tmpvisu.     , 2, 2),
                                                    CvColor.Red);
                                }

                            }


                        } 
                        tmpvisu . SaveImage ( „10-visu.png” ); 
                        tmpvisu . Zapisać obraz("10-visu" + file + ".png");
                    }
                    #endregion

                }

            } Konsola . WriteLine (
            "total error: {0}, ideal Pixel per Grain: {1:0.0}", totalError, totalPixels*1.0/expectedGrains.Sum());

        }
    }
}

Mały test wykorzystujący zakodowany na stałe rozmiar pikseli na ziarno wynoszący 2544,4 wykazał całkowity błąd 36, który jest nadal większy niż większość innych rozwiązań.

wprowadź opis zdjęcia tutaj wprowadź opis zdjęcia tutaj wprowadź opis zdjęcia tutaj wprowadź opis zdjęcia tutaj

HugoRune
źródło
Myślę, że można użyć progu (operacja erodowania może być również przydatna) z pewną niewielką wartością wyniku transformacji odległości - to powinno podzielić niektóre grupy ziaren na mniejsze grupy (najlepiej - tylko z 1 lub 2 ziarnami). Niż powinno być łatwiej policzyć te samotne ziarna. Duże grupy, które można zaliczyć jako większość ludzi tutaj - dzieląc powierzchnię przez średnią powierzchnię pojedynczego ziarna.
cyriel,
9

HTML + JavaScript: wynik 39

Dokładne wartości to:

Estimated | Actual
        3 |      3
        5 |      5
       12 |     12
       23 |     25
       51 |     50
       82 |     83
      125 |    120
      161 |    150
      167 |    151
      223 |    200

Rozkłada się (nie jest dokładny) na większe wartości.

window.onload = function() {
  var $ = document.querySelector.bind(document);
  var canvas = $("canvas"),
    ctx = canvas.getContext("2d");

  function handleFileSelect(evt) {
    evt.preventDefault();
    var file = evt.target.files[0],
      reader = new FileReader();
    if (!file) return;
    reader.onload = function(e) {
      var img = new Image();
      img.onload = function() {
        canvas.width = this.width;
        canvas.height = this.height;
        ctx.drawImage(this, 0, 0);
        start();
      };
      img.src = e.target.result;
    };
    reader.readAsDataURL(file);
  }


  function start() {
    var imgdata = ctx.getImageData(0, 0, canvas.width, canvas.height);
    var data = imgdata.data;
    var background = 0;
    var totalPixels = data.length / 4;
    for (var i = 0; i < data.length; i += 4) {
      var red = data[i],
        green = data[i + 1],
        blue = data[i + 2];
      if (Math.abs(red - 197) < 40 && Math.abs(green - 176) < 40 && Math.abs(blue - 133) < 40) {
        ++background;
        data[i] = 1;
        data[i + 1] = 1;
        data[i + 2] = 1;
      }
    }
    ctx.putImageData(imgdata, 0, 0);
    console.log("Pixels of rice", (totalPixels - background));
    // console.log("Total pixels", totalPixels);
    $("output").innerHTML = "Approximately " + Math.round((totalPixels - background) / 2670) + " grains of rice.";
  }

  $("input").onchange = handleFileSelect;
}
<input type="file" id="f" />
<canvas></canvas>
<output></output>

Objaśnienie: Zasadniczo zlicza liczbę pikseli ryżu i dzieli ją przez średnią liczbę pikseli na ziarno.

soktinpk
źródło
Używając obrazu z 3
ryżami
1
@Kroltan Nie, gdy używasz obrazu w pełnym rozmiarze .
Calvin's Hobbies
1
@ Calvin'sHobbies FF36 na Windows dostaje 0, na Ubuntu dostaje 3, z obrazem w pełnym rozmiarze.
Kroltan
4
@BobbyJack Ryż ma gwarantowaną mniej więcej taką samą skalę na wszystkich obrazach. Nie widzę z tym żadnych problemów.
Calvin's Hobbies
1
@githubphagocyte - wyjaśnienie jest dość oczywiste - jeśli policzysz wszystkie białe piksele w wyniku binaryzacji obrazu i podzielisz tę liczbę przez liczbę ziaren w obrazie, otrzymasz ten wynik. Oczywiście dokładny wynik może się różnić, ze względu na zastosowaną metodę binaryzacji i inne rzeczy (np. Operacje wykonywane po binaryzacji), ale jak widać w innych odpowiedziach, będzie ona w zakresie 2500-3500.
cyriel,
4

Próba z php. Nie odpowiedź o najniższej punktacji, ale dość prosty kod

WYNIK: 31

<?php
for($c = 1; $c <= 10; $c++) {
  $a = imagecreatefromjpeg("/tmp/$c.jpg");
  list($width, $height) = getimagesize("/tmp/$c.jpg");
  $rice = 0;
  for($i = 0; $i < $width; $i++) {
    for($j = 0; $j < $height; $j++) {
      $colour = imagecolorat($a, $i, $j);
      if (($colour & 0xFF) < 95) $rice++;
    }
  }
  echo ceil($rice/2966);
}

Samoocena

95 to niebieska wartość, która wydawała się działać, gdy testowanie za pomocą GIMP 2966 jest średnią wielkością ziarna

exussum
źródło