User:Swedmann/R
From RationalWiki
#!/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)