User:Swedmann/R

#!/usr/bin/Rscript
require (biOps)
require (RPostgreSQL)
require (RODBC)

dsource = 'postgres'
dbname = 'guestdb'
duid = 'guest'
dpwd = 'guest'

pg = PostgreSQL ()
pgconn = dbConnect (pg, user = duid, password = dpwd, dbname = dbname)
conn = odbcConnect (dsource, uid = duid, pwd = dpwd)

gw_fft = list ()
for (n in 1:40) { gw_fft[[n]] = c (0); dim (gw_fft[[n]]) = 1 }

readimg = function (file)
  {
    img = readJpeg (file)

    if (length (dim (img)) == 3)
      {
        img = img[,,1] + img[,,2] + img[,,3]
        img = imagedata (img)
      }

    imgNormalize (img)
  }

makegw = function (pgconn, conn, d)
  {
    if (all (dim (gw_fft[[1]]) == d)) return (gw_fft)

    gw = dbGetQuery (pgconn, paste ('select gw from gw where width = ',
      d[2], ' and height = ', d[1], ';', sep = ''))

    if (length (gw) != 0)
      {
        gw_fft = unserialize (charToRaw (gw$gw[1]))
        return (gw_fft)
      }

    gw_fft <<- list ()
    s = .5
    s2 = s ** 2
    n = 1

    for (rot in pi / 8 * 0:7)
      {
        for (scl in 2 ** (1:5))
          {
            gw = array (0, d)
            k = scl * c (cos (rot), sin (rot))
            k2 = sum (k * k)

            for (j in 1:d[1])
              {
                for (i in 1:d[2])
                  {
                    J = j / d[1] - .5
                    I = i / d[2] - .5
                    x = c (I, J)
                    x2 = sum (x * x)

                    gw[j, i] = k2 / s2 * exp (- k2 * x2 / (2 * s2)) *
                      (exp (1i * sum (k * x)) - exp (-s2 / 2))
                  }
              }

            gw_fft[[n]] <<- imagedata (imgFFT (imagedata (imgFFTShift (gw))))
            n = n + 1
          }
      }

    entry = data.frame (width = d[2], height = d[1],
      gw = rawToChar (serialize (gw_fft, NULL, ascii = T)))
    sqlSave (conn, entry, tablename = 'gw', rownames = F, append = T, fast = F)

    gw_fft
  }

getfeat = function (pgconn, conn, imgdata)
  {
    gwfft = makegw (pgconn, conn, dim (imgdata))
    fft = imagedata (imgFFT (imgdata))
    a = list ()
    w = 9

    for (n in 1:length (gwfft))
      a[[n]] = imgFFTInv (fft * gwfft[[n]] * 2 ** 10) / 2 ** 10

    mask = array (0, dim (a[[1]]))

    for (n in 1:length (a))
      {
        mask[(a[[n]] == imgMaximumFilter (a[[n]], w)) & (a[[n]]
                > imgConvolve (a[[n]], array (1 / w ** 2, c (w, w))))] = 1
      }

    feat = c ()

    ya = as.integer (w / 2)
    yb = as.integer (dim (mask)[1] - w / 2)
    xa = as.integer (w / 2)
    xb = as.integer (dim (mask)[2] - w / 2)

    for (j in ya:yb)
      {
        for (i in xa:xb)
          {
            if (mask[j, i])
              {
                feat = append (feat, c (i / dim (mask)[2],
                  j / dim (mask)[1]))

                for (n in 1:length (a))
                  feat = append (feat, a[[n]][j, i])
              }
          }
      }

    feat
  }

featcmp = function (a, b, dist, thold)
  {
    count = 0
    sim = 0

    for (i in 42 * 1:(length (a) / 42) - 41)
      {
        s = 0

        for (j in 42 * 1:(length (b) / 42) - 41)
          {
            d = a[i:(i + 1)] - b[j:(j + 1)]
            d = sum (d * d)

            if (d < dist * dist)
              {
                nom = sum (a[(i + 2):(i + 41)] * b[(j + 2):(j + 41)])
                denom = sqrt (sum (a[(i + 2):(i + 41)] ** 2)
                  * sum (b[(j + 2):(j + 41)] ** 2));

                if (denom != 0 && nom / denom > s)
                  s = nom / denom
              }
          }

        if (s > thold)
          {
            sim = sim + s
            count = count + 1
          }
      }

    c (sim, count, length (a) / 42, length (b) / 42)
  }

add_new = function (pgconn, conn, imglist, name, data)
  {
    feat = c ()

    for (i in imglist)
      {
        feat = append (feat, getfeat (pgconn, conn, i))
      }

    entry = data.frame (feat = rawToChar (serialize (feat, NULL, ascii = T)),
      data = data, name = name)

    sqlSave (conn, entry, tablename = 'frfr', rownames = F,
             append = T, fast = F)
  }

match = function (pgconn, conn, img, dist, thold)
  {
    feat = getfeat (pgconn, conn, img)
    fdata = dbReadTable (pgconn, 'frfr', row.names = F)
    sim = c ()

    for (i in 1:length (fdata$feat))
      {
        s = featcmp (feat, unserialize (charToRaw (fdata$feat[i])),
          dist, thold)
        sim[i] = s[1] / s[2]
      }

    o = order (sim, decreasing = T)
    data.frame (name = fdata$name[o], data = fdata$data[o], sim = sim[o])
  }

buildfeat = function (pgconn, conn, img, w)
  {
    feat = list ()

    for (i in 1:length (img))
      {
        feat[[i]] = getfeat (pgconn, conn, img[[i]][[w[1]]])

        for (j in w[-1])
          {
            feat[[i]] = append (feat[[i]], getfeat (pgconn, conn,
                  img[[i]][[j]]))
          }
      }

    feat
  }

featmat = function (a, b, dist, thold)
  {
    s = array (0, c (length (a), length (b), 4))

    for (i in 1:length (a))
      {
        for (j in 1:length (b))
          {
            sim = featcmp (a[[i]], b[[j]], dist, thold)
            s[i, j, ] = sim
          }
      }

    s
  }

arg = commandArgs (trailingOnly = T)

if (arg[1] == 'add')
  {
    name = arg[2]
    data = arg[3]
    imgfile = arg[4:length (arg)]

    img = list ()

    for (i in 1:length (imgfile))
      {
        img[[i]] = readimg (imgfile[i])
      }

    add_new (pgconn, conn, img, name, data)
  }

if (arg[1] == 'match')
  {
    dist = as.double (arg[2])
    thold = as.double (arg[3])
    img = readimg (arg[4])

    s = match (pgconn, conn, img, dist, thold)

    write (s$name[1], file = stdout ())
    write (s$data[1], file = stdout ())
    write (s$sim[1], file = stdout ())
  }

odbcClose (conn)
dbDisconnect (pgconn)