Contando granos de arroz

81

Considere estas 10 imágenes de varias cantidades de granos crudos de arroz blanco.
ESTAS SON SOLO MINIATURAS. Haga clic en una imagen para verla a tamaño completo.

A: B: C: D: E:UNA si C re mi

F: G: H: I: J:F sol H yo J

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

Darse cuenta de...

  • Los granos pueden tocarse entre sí, pero nunca se superponen. El diseño de los granos nunca tiene más de un grano de alto.
  • Las imágenes tienen diferentes dimensiones, pero la escala del arroz en todas ellas es consistente porque la cámara y el fondo eran estacionarios.
  • Los granos nunca salen de los límites ni tocan los límites de la imagen.
  • El fondo es siempre el mismo tono consistente de blanco amarillento.
  • Los granos pequeños y grandes se cuentan por igual como un grano cada uno.

Estos 5 puntos son garantías para todas las imágenes de este tipo.

Desafío

Escriba un programa que tome esas imágenes y, con la mayor precisión posible, cuente la cantidad de granos de arroz.

Su programa debe tomar el nombre de archivo de la imagen e imprimir la cantidad de granos que calcula. Su programa debe funcionar al menos para uno de estos formatos de archivo de imagen: JPEG, mapa de bits, PNG, GIF, TIFF (en este momento las imágenes son todas JPEG).

Usted puede utilizar las bibliotecas de procesamiento de imágenes y visión por ordenador.

No puede codificar las salidas de las 10 imágenes de ejemplo. Su algoritmo debe ser aplicable a todas las imágenes similares de grano de arroz. Debería poder ejecutarse en menos de 5 minutos en una computadora moderna decente si el área de la imagen es inferior a 2000 * 2000 píxeles y hay menos de 300 granos de arroz.

Puntuación

Para cada una de las 10 imágenes, tome el valor absoluto de la cantidad real de granos menos la cantidad de granos que predice su programa. Suma estos valores absolutos para obtener tu puntaje. El puntaje más bajo gana. Un puntaje de 0 es perfecto.

En caso de empate, gana la respuesta más votada. Puedo probar su programa en imágenes adicionales para verificar su validez y precisión.

Pasatiempos de Calvin
fuente
1
¡Seguramente alguien tiene que probar scikit-learn!
Gran concurso! :) Por cierto, ¿podría contarnos algo sobre la fecha de finalización de este desafío?
cyriel
1
@Lembik Hasta 7 :)
Dr. belisarius
55
Un día, un científico de arroz vendrá y se pondrá loco de felicidad por la existencia de esta pregunta.
Nit
2
@Nit Solo diles ncbi.nlm.nih.gov/pmc/articles/PMC3510117 :)
Dr. belisarius

Respuestas:

22

Mathematica, puntuación: 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;

Creo que los nombres de las funciones son lo suficientemente descriptivos:

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]

Procesando todas las imágenes a la vez:

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

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

El puntaje es:

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

Aquí puede ver la sensibilidad de la puntuación con el tamaño de grano utilizado:

Gráficos de Mathematica

Dr. belisario
fuente
2
Mucho más claro, gracias!
Aficiones de Calvin
¿Se puede copiar este procedimiento exacto en Python o hay algo especial que Mathematica está haciendo aquí que las bibliotecas de Python no pueden hacer?
@Lembik No tengo idea. No hay pitón aquí. Lo siento. (Sin embargo, dudo exactamente los mismos algoritmos para EdgeDetect[], DeleteSmallComponents[]y Dilation[]se implementan en otro lugar)
Dr. belisarius
55

Python, Puntuación: 24 dieciséis

Esta solución, como la de Falko, se basa en medir el área de "primer plano" y dividirla por el área de grano promedio.

De hecho, lo que este programa intenta detectar es el fondo, no tanto como el primer plano. Usando el hecho de que los granos de arroz nunca tocan el límite de la imagen, el programa comienza rellenando con blanco en la esquina superior izquierda. El algoritmo de relleno de inundación pinta píxeles adyacentes si la diferencia entre ellos y el brillo del píxel actual se encuentra dentro de un cierto umbral, ajustándose así a un cambio gradual en el color de fondo. Al final de esta etapa, la imagen podría verse así:

Figura 1

Como puede ver, hace un buen trabajo al detectar el fondo, pero deja fuera cualquier área que esté "atrapada" entre los granos. Manejamos estas áreas al estimar el brillo del fondo en cada píxel y al paitizar todos los píxeles iguales o más brillantes. Esta estimación funciona así: durante la etapa de inundación, calculamos el brillo de fondo promedio para cada fila y cada columna. El brillo de fondo estimado en cada píxel es el promedio del brillo de la fila y la columna en ese píxel. Esto produce algo como esto:

Figura 2

EDITAR: Finalmente, el área de cada región de primer plano continuo (es decir, no blanca) se divide por el área de grano promedio, precalculada, lo que nos da una estimación del recuento de granos en dicha región. La suma de estas cantidades es el resultado. Inicialmente, hicimos lo mismo para toda el área de primer plano en su conjunto, pero este enfoque es, literalmente, más detallado.


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

Toma el nombre de archivo de entrada a través de la línea de comando.

Resultados

      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

UNA si C re mi

F sol H yo J

Ana
fuente
2
Esta es una solución realmente inteligente, ¡buen trabajo!
Chris Cirefice
1
de donde avg_grain_area = 3038.38;viene
njzk2
2
¿eso no cuenta como hardcoding the result?
njzk2
55
@ njzk2 No. Dada la regla The images have different dimensions but the scale of the rice in all of them is consistent because the camera and background were stationary.Esto es simplemente un valor que representa esa regla. El resultado, sin embargo, cambia según la entrada. Si cambia la regla, este valor cambiará, pero el resultado será el mismo, según la entrada.
Adam Davis el
66
Estoy bien con lo del área promedio. El área de grano es (aproximadamente) constante en las imágenes.
Aficiones de Calvin
28

Python + OpenCV: puntaje 27

Escaneo de linea horizontal

Idea: escanee la imagen, una fila a la vez. Para cada línea, cuente el número de granos de arroz encontrados (verificando si el píxel se vuelve negro a blanco o lo contrario). Si el número de granos para la línea aumenta (en comparación con la línea anterior), significa que encontramos un nuevo grano. Si ese número disminuye, significa que pasamos por alto un grano. En este caso, agregue +1 al resultado total.

ingrese la descripción de la imagen aquí

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

Debido a la forma en que funciona el algoritmo, es importante tener una imagen limpia en blanco y negro. Mucho ruido produce malos resultados. Primero se limpia el fondo principal usando el relleno de inundación (solución similar a la respuesta de Ell), luego se aplica el umbral para producir un resultado en blanco y negro.

ingrese la descripción de la imagen aquí

Está lejos de ser perfecto, pero produce buenos resultados con respecto a la simplicidad. Probablemente haya muchas maneras de mejorarlo (proporcionando una mejor imagen en blanco y negro, escaneando en otras direcciones (por ejemplo: vertical, diagonal) tomando el promedio, etc.)

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

Los errores por imagen: 0, 0, 0, 3, 0, 12, 4, 0, 7, 1

tigrou
fuente
24

Python + OpenCV: Puntuación 84

Aquí hay un primer intento ingenuo. Aplica un umbral adaptativo con parámetros ajustados manualmente, cierra algunos agujeros con posterior erosión y dilución y deriva el número de granos del área de primer plano.

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)

Aquí puede ver las imágenes binarias intermedias (el negro es primer plano):

ingrese la descripción de la imagen aquí

Los errores por imagen son 0, 0, 2, 2, 4, 0, 27, 42, 0 y 7 granos.

Falko
fuente
20

C # + OpenCvSharp, Puntuación: 2

Este es mi segundo intento. Es bastante diferente de mi primer intento , que es mucho más simple, por lo que lo publico como una solución separada.

La idea básica es identificar y etiquetar cada grano individual mediante un ajuste de elipse iterativo. Luego, elimine los píxeles de este grano de la fuente e intente encontrar el siguiente grano, hasta que se haya etiquetado cada píxel.

Esta no es la solución más bonita. Es un cerdo gigante con 600 líneas de código. Necesita 1,5 minutos para la imagen más grande. Y realmente me disculpo por el código desordenado.

Hay tantos parámetros y formas de pensar en esto que tengo mucho miedo de sobreajustar mi programa para las 10 imágenes de muestra. El puntaje final de 2 es casi definitivamente un caso de sobreajuste: tengo dos parámetros, average grain size in pixely minimum ratio of pixel / elipse_area, al final, simplemente agoté todas las combinaciones de estos dos parámetros hasta que obtuve el puntaje más bajo. No estoy seguro de si esto es tan kosher con las reglas de este desafío.

average_grain_size_in_pixel = 2530
pixel / elipse_area >= 0.73

Pero incluso sin estos embragues sobreajustados, los resultados son bastante buenos. Sin un tamaño de grano fijo o una relación de píxeles, simplemente estimando el tamaño de grano promedio a partir de las imágenes de entrenamiento, el puntaje sigue siendo 27.

Y obtengo como salida no solo el número, sino la posición, orientación y forma reales de cada grano. hay una pequeña cantidad de granos mal etiquetados, pero en general la mayoría de las etiquetas coinciden con precisión con los granos reales:

A UNA B si C C D re Emi

F F G sol H H I yo JJ

(haga clic en cada imagen para la versión de tamaño completo)

Después de este paso de etiquetado, mi programa analiza cada grano individual y las estimaciones basadas en el número de píxeles y la relación píxel / área de elipse, ya sea

  • un solo grano (+1)
  • granos múltiples mal etiquetados como uno (+ X)
  • una gota demasiado pequeña para ser un grano (+0)

Los puntajes de error para cada imagen son A:0; B:0; C:0; D:0; E:2; F:0; G:0 ; H:0; I:0, J:0

Sin embargo, el error real es probablemente un poco más alto. Algunos errores dentro de la misma imagen se cancelan mutuamente. La imagen H en particular tiene algunos granos mal etiquetados, mientras que en la imagen E las etiquetas son en su mayoría correctas

El concepto es un poco artificial:

  • Primero, el primer plano se separa a través del umbral de otsu en el canal de saturación (vea mi respuesta anterior para más detalles)

  • repita hasta que no queden más píxeles:

    • seleccione la gota más grande
    • elija 10 píxeles de borde aleatorio en este blob como posiciones iniciales para un grano

    • para cada punto de partida

      • Asumir un grano con altura y ancho de 10 píxeles en esta posición.

      • repetir hasta la convergencia

        • vaya radialmente hacia afuera desde este punto, en diferentes ángulos, hasta que encuentre un píxel de borde (blanco a negro)

        • Con suerte, los píxeles encontrados deberían ser los píxeles de borde de un solo grano. Intente separar los datos internos de los valores atípicos, descartando píxeles que estén más alejados de la elipse supuesta que los demás.

        • trate repetidamente de ajustar una elipse a través de un subconjunto de los signos internos, mantenga la mejor elipse (RANSACK)

        • Actualice la posición del grano, la orientación, el ancho y la altura con la elipse encontrada

        • Si la posición del grano no cambia significativamente, deténgase

    • entre los 10 granos ajustados, elija el mejor grano de acuerdo con la forma, el número de píxeles de borde. Descarta a los demás

    • elimine todos los píxeles de este grano de la imagen de origen, luego repita

    • finalmente, revise la lista de granos encontrados y cuente cada grano como 1 grano, 0 granos (demasiado pequeño) o 2 granos (demasiado grande)

Uno de mis principales problemas fue que no quería implementar una métrica de distancia de punto de elipse completa, ya que calcular eso en sí mismo es un proceso iterativo complicado. Así que utilicé varias soluciones usando las funciones de OpenCV Ellipse2Poly y FitEllipse, y los resultados no son demasiado bonitos.

Aparentemente también rompí el límite de tamaño para codegolf.

Una respuesta está limitada a 30000 caracteres, actualmente estoy en 34000. Así que tendré que acortar un poco el código a continuación.

El código completo se puede ver en http://pastebin.com/RgM7hMxq

Lo siento, no sabía que había un límite de tamaño.

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);
    }

}

Estoy un poco avergonzado con esta solución porque a) no estoy seguro de si está dentro del espíritu de este desafío, yb) es demasiado grande para una respuesta de codegolf y carece de la elegancia de las otras soluciones.

Por otro lado, estoy bastante contento con el progreso que logré en el etiquetado de los granos, no solo en contarlos, así que eso es todo.

HugoRune
fuente
Usted sabe que puede reducir la longitud de ese código en magnitudes usando nombres más pequeños y aplicando algunas otras técnicas de golf;)
Optimizer
Probablemente, pero no quería ofuscar aún más esta solución. Es demasiado ofuscado para mis gustos, ya que es :)
HugoRune
+1 por el esfuerzo y porque eres el único que encuentra la manera de mostrar individualmente cada grano. Desafortunadamente, el código está un poco hinchado y depende mucho de las constantes codificadas. Me gustaría ver cómo funciona el algoritmo de línea de exploración que escribí en esto (en los granos de color inviduales).
tigrou
Realmente creo que este es el enfoque correcto para este tipo de problema (+1 para usted), pero una cosa que me pregunto es por qué "elige 10 píxeles de borde aleatorio", creo que obtendrá un mejor rendimiento si elige los puntos de borde con el menor número de puntos de borde cercanos (es decir, partes que sobresalen), creo que (en teoría) esto eliminaría primero los granos "más fáciles", ¿ha considerado esto?
David Rogers
Lo he pensado, pero aún no lo he probado. La "posición inicial aleatoria 10" fue una adición tardía, que fue fácil de agregar y de paralelizar. Antes de eso, 'una posición inicial aleatoria' era mucho mejor que 'siempre la esquina superior izquierda'. El peligro de elegir las posiciones iniciales con la misma estrategia cada vez es que cuando elimine el mejor ajuste, los otros 9 probablemente serán elegidos nuevamente la próxima vez, y con el tiempo la peor de esas posiciones iniciales se quedará atrás y será elegida nuevamente y de nuevo. Una parte que sobresale puede ser solo los restos de un grano anterior eliminado de forma incompleta.
HugoRune
17

C ++, OpenCV, puntaje: 9

La idea básica de mi método es bastante simple: intente borrar granos individuales (y "granos dobles" - 2 (¡pero no más!) Granos, cerca uno del otro) de la imagen y luego cuente el descanso usando el método basado en el área (como Falko, Ell y Belisario). Usar este enfoque es un poco mejor que el "método de área" estándar, porque es más fácil encontrar un buen valor promedio de PixelesPerObjeto.

(1er paso) Primero que nada necesitamos usar la binarización de Otsu en el canal S de imagen en HSV. El siguiente paso es usar el operador de dilatación para mejorar la calidad del primer plano extraído. De lo que necesitamos encontrar contornos. Por supuesto, algunos contornos no son granos de arroz: necesitamos eliminar contornos demasiado pequeños (con un área más pequeña que el promedio de píxeles por objeto / 4. El promedio de píxeles por objeto es 2855 en mi situación). Ahora finalmente podemos comenzar a contar granos :) (2do paso) Encontrar granos simples y dobles es bastante simple: solo busque en la lista de contornos los contornos con área dentro de rangos específicos; si el área de contorno está dentro del rango, elimínelo de la lista y agregue 1 (o 2 si era "doble" de grano) al contador de granos. (3er paso) El último paso es, por supuesto, dividir el área de los contornos restantes por el valor promedio de PixelPerObject y agregar el resultado al contador de granos.

Las imágenes (para la imagen F.jpg) deberían mostrar esta idea mejor que las palabras:
primer paso (sin contornos pequeños (ruido)): 1er paso (sin contornos pequeños (ruido))
segundo paso - solo contornos simples: 2do paso - solo contornos simples
tercer paso - contornos restantes: 3er paso - contornos restantes

Aquí está el código, es bastante feo, pero debería funcionar sin ningún problema. Por supuesto, se requiere 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;
}

Si desea ver los resultados de todos los pasos, elimine el comentario de todas las llamadas a funciones imshow (.., ..) y establezca la variable fastProcessing en false. Las imágenes (A.jpg, B.jpg, ...) deben estar en las imágenes del directorio. Alternativamente, puede dar el nombre de una imagen como parámetro desde la línea de comandos.

Por supuesto, si algo no está claro, puedo explicarlo y / o proporcionar algunas imágenes / información.

cyriel
fuente
12

C # + OpenCvSharp, puntuación: 71

Esto es muy irritante, traté de obtener una solución que realmente identifique cada grano usando la cuenca hidrográfica , pero simplemente. hipocresía. obtener. eso. a. trabajo.

Me decidí por una solución que al menos separa algunos granos individuales y luego usa esos granos para estimar el tamaño promedio de grano. Sin embargo, hasta ahora no puedo superar las soluciones con un tamaño de grano codificado.

Entonces, lo más destacado de esta solución: no supone un tamaño de píxel fijo para los granos, y debería funcionar incluso si se mueve la cámara o se cambia el tipo de arroz.

A.jpg; cantidad de granos: 3; esperado 3; error 0; píxeles por grano: 2525,0;
B.jpg; cantidad de granos: 7; esperado 5; error 2; píxeles por grano: 1920,0;
C.jpg; cantidad de granos: 6; esperado 12; error 6; píxeles por grano: 4242,5;
D.jpg; cantidad de granos: 23; esperado 25; error 2; píxeles por grano: 2415,5;
E.jpg; cantidad de granos: 47; esperado 50; error 3; píxeles por grano: 2729,9;
F.jpg; cantidad de granos: 65; esperado 83; error 18; píxeles por grano: 2860,5;
G.jpg; cantidad de granos: 120; esperado 120; error 0; píxeles por grano: 2552,3;
H.jpg; cantidad de granos: 159; esperado 150; error 9; píxeles por grano: 2624,7;
I.jpg; cantidad de granos: 141; esperado 151; error 10; píxeles por grano: 2697,4;
J.jpg; cantidad de granos: 179; esperado 200; error 21; píxeles por grano: 2847,1;
error total: 71

Mi solución funciona así:

Separe el primer plano transformando la imagen a HSV y aplicando el umbral de Otsu en el canal de saturación. Esto es muy simple, funciona extremadamente bien, y lo recomendaría a todos los que quieran probar este desafío:

saturation channel                -->         Otsu thresholding

ingrese la descripción de la imagen aquí -> ingrese la descripción de la imagen aquí

Esto eliminará limpiamente el fondo.

Luego, además, eliminé las sombras de grano del primer plano, aplicando un umbral fijo al canal de valores. (No estoy seguro si eso realmente ayuda mucho, pero fue lo suficientemente simple como para agregar)

ingrese la descripción de la imagen aquí

Luego aplico una transformación de distancia en la imagen de primer plano.

ingrese la descripción de la imagen aquí

y encuentre todos los máximos locales en esta transformación de distancia.

Aquí es donde se rompe mi idea. para evitar obtener máximos locales múltiples dentro del mismo grano, tengo que filtrar mucho. Actualmente mantengo solo el máximo más fuerte dentro de un radio de 45 píxeles, lo que significa que no todos los granos tienen un máximo local. Y realmente no tengo una justificación para el radio de 45 píxeles, solo fue un valor que funcionó.

ingrese la descripción de la imagen aquí

(Como puede ver, esas no son semillas suficientes para dar cuenta de cada grano)

Luego uso esos máximos como semillas para el algoritmo de cuenca:

ingrese la descripción de la imagen aquí

Los resultados son meh . Esperaba principalmente granos individuales, pero los grupos todavía son demasiado grandes.

Ahora identifico los blobs más pequeños, cuento su tamaño de píxel promedio y luego calculo el número de granos a partir de eso. Esto no es lo que planeé hacer al principio, pero esta era la única forma de salvarlo.

utilizando el sistema ; 
utilizando el sistema . Colecciones . Genérico ; 
utilizando el sistema . Linq ; 
utilizando el sistema . Texto ; 
usando OpenCvSharp ;

espacio de nombres GrainTest2 { Programa de clase { vacío estático Principal ( cadena [] args ) { cadena [] archivos = nuevo [] { "sourceA.jpg" , "sourceB.jpg" , "sourceC.jpg" , "sourceD.jpg" , " sourceE.jpg " , " sourceF.jpg " , " sourceG.jpg " , " sourceH.jpg " , " sourceI.jpg " , " sourceJ.jpg " , };int [] esperadaGrains

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

            int totalError = 0 ; int totalPixels = 0 ; 
             

            for ( int fileNo = 0 ; fileNo markers = new List (); 
                    using ( CvMemStorage storage = new CvMemStorage ()) 
                    using ( CvContourScanner scanner = new CvContourScanner ( localMaxima , storage , CvContour . SizeOf , ContourRetrieval . External , ContourChain) . ApproxNone ).         
                    { // establece cada máximo local como semilla número 25, 35, 45, ... // (los números reales no importan, elegidos para una mejor visibilidad en png) int markerNo = 20 ; foreach ( CvSeq c en el escáner ) { 
                            markerNo + = 5 ; 
                            marcadores . Añadir ( marcador No ); 
                            waterShedMarkers . DrawContours ( c , nuevo CvScalar ( marcador No ), nuevo
                        
                        
                         
                         
                             CvScalar ( marcador No ), 0 , - 1 ); } } 
                    waterShedMarkers . SaveImage ( "08-watershed-seeds.png" );  
                        
                    


                    fuente . Watershed ( waterShedMarkers ); 
                    waterShedMarkers . SaveImage ( "09-watershed-result.png" );


                    Lista píxelesPerBlob = nueva Lista ();  

                    // Terrible hack because I could not get Cv2.ConnectedComponents to work with this openCv wrapper
                    // So I made a workaround to count the number of pixels per blob
                    waterShedMarkers.ConvertScale(waterShedThreshold);
                    foreach (int markerNo in markers)
                    {
                        using (CvMat tmp = new CvMat(waterShedMarkers.Rows, waterShedThreshold.Cols, MatrixType.U8C1))
                        {
                            waterShedMarkers.CmpS(markerNo, tmp, ArrComparison.EQ);
                            pixelsPerBlob.Add(tmp.CountNonZero());

                        }
                    }

                    // estimate the size of a single grain
                    // step 1: assume that the 10% smallest blob is a whole grain;
                    double 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.External,
                                                                                ContourChain.ApproxNone))
                            {
                                tmpvisu.Set(CvColor.Random(), tmp);
                                foreach (CvSeq c in scanner)
                                {
                                    //tmpvisu.DrawContours(c, CvColor.Random(), CvColor.DarkGreen, 0, -1);
                                    tmpvisu.PutText("" + Math.Round(curGrains, 1), c.First().Value, new CvFont(FontFace.HersheyPlain, 2, 2),
                                                    CvColor.Red);
                                }

                            }


                        }
                        tmpvisu.SaveImage("10-visu.png");
                        tmpvisu.SaveImage("10-visu" + file + ".png");
                    }
                    #endregion

                }

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

        }
    }
}

Una pequeña prueba con un tamaño de píxel por grano codificado de 2544.4 mostró un error total de 36, que aún es más grande que la mayoría de las otras soluciones.

ingrese la descripción de la imagen aquí ingrese la descripción de la imagen aquí ingrese la descripción de la imagen aquí ingrese la descripción de la imagen aquí

HugoRune
fuente
Creo que puede usar el umbral (la operación de erosión también podría ser útil) con algún valor pequeño en el resultado de la transformación de distancia; esto debería dividir algunos grupos de granos en grupos más pequeños (preferiblemente, con solo 1 o 2 granos). De lo que debería ser más fácil contar esos granos solitarios. Grandes grupos que puede contar como la mayoría de las personas aquí, dividiendo el área por el área promedio de grano único.
cyriel
9

HTML + Javascript: Puntuación 39

Los valores exactos son:

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

Se descompone (no es preciso) en los valores más grandes.

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>

Explicación: Básicamente, cuenta el número de píxeles de arroz y lo divide por el promedio de píxeles por grano.

soktinpk
fuente
Usando la imagen de 3 arroces, estimó 0 para mí ...: /
Kroltan
1
@Kroltan No cuando usa la imagen a tamaño completo .
Las aficiones de Calvin el
1
@ Calvin'sHobbies FF36 en Windows obtiene 0, en Ubuntu obtiene 3, con la imagen a tamaño completo.
Kroltan el
44
@BobbyJack Se garantiza que el arroz tendrá más o menos la misma escala en las imágenes. No veo problemas con eso.
Aficiones de Calvin
1
@githubphagocyte: una explicación es bastante obvia: si cuenta todos los píxeles blancos en el resultado de la binarización de la imagen y divide este número por el número de granos en la imagen, obtendrá este resultado. Por supuesto, el resultado exacto puede diferir, debido al método de binarización utilizado y otras cosas (como las operaciones realizadas después de la binarización), pero como puede ver en otras respuestas, estará en el rango 2500-3500.
cyriel
4

Un intento con php, no la respuesta con la puntuación más baja, sino su código bastante simple

PUNTUACIÓN: 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);
}

Auto puntuación

95 es un valor azul que parece funcionar cuando se prueba con GIMP 2966 es el tamaño de grano promedio

exussum
fuente