Kristian/scanned-image-extractor
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
v0.2.601-beta
scanned-image-extractor/scannerExtract/preloadsource.cpp

1966 lines
67 KiB
C++
Raw Permalink Normal View History

initial commit
2024-01-29 23:24:37 +00:00
/***********************************************************************
* This file is part of Scanned Image Extract.
*
* Scanned Image Extract is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Scanned Image Extract is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Scanned Image Extract. If not, see <http://www.gnu.org/licenses/>
*
*
* Copyright (C) 2015, Dominik Rueß; info@dominik-ruess.de
**********************************************************************/
#include "preloadsource.h"
#include <lbfgs.h>
#include <QMutexLocker>
#include <QDebug>
#include <QtMath>
#include <QDateTime>
#include <opencv/cv.h>
#ifndef OPENCV2
#include <opencv2/opencv.hpp>
#include <opencv2/imgcodecs.hpp>
#include <opencv2/imgproc.hpp>
#else
#include <opencv/highgui.h>
#endif
#include <vector>
//#define DEBUG
//#define DEBUG_TEST
#include "copytargets.h"
#define CvPt2QtPt(p) QPointF(p.x, p.y)
QMutex PreloadSource::preloadMutex;
QWaitCondition PreloadSource::waitFinished;
#ifdef DEBUG_TEST
static int tmpVal2 = 1;
#endif
#ifdef DEBUG
static int tmpVal =0;
#endif
#define ABORT_COMPUTATION(s) \
if (_abortAndLoadNext || _stopped) return s;
PreloadSource::PreloadSource(QObject *parent)
: QThread(parent)
, _stopped(false)
, _abortAndLoadNext(false)
, _qImageFmt(QImage::Format_ARGB32)
{
//setEnergyWeights(5,1,1,5);
}
PreloadSource::~PreloadSource()
{
_fileListMutex.lock();
_filesToLoad.clear();
_stopped = true;
_condition.wakeOne();
_fileListMutex.unlock();
wait();
}
void PreloadSource::addLoadFiles(const SourceFilePtr& file,
const int position,
const bool highPriority,
const bool onlyLoadImages)
{
QMutexLocker l(&_fileListMutex);
if (_current && _current->id == file->id) {
return;
}
if (highPriority) {
_abortAndLoadNext = true;
}
_filesToLoad.push_front(QPair<SourceFilePtr, int> (file, onlyLoadImages ? -position-1 : position)); // negative position: onlyLoadImages
// now remove duplicates (from back)
for (int i=_filesToLoad.size()-1; i>= 0 ; i--) {
for (int j=i-1; j >=0; j--) {
if (_filesToLoad[i].first->source.fileName() == _filesToLoad[j].first->source.fileName()) {
_filesToLoad.erase(_filesToLoad.begin() + i);
break;
}
}
}
_stopped = false;
if (!isRunning()) {
if (highPriority) {
start(QThread::HighPriority);
}
else
{
start(QThread::LowPriority);
}
} else {
if (highPriority) {
setPriority(QThread::HighPriority);
}
else
{
setPriority(QThread::LowPriority);
}
_condition.wakeOne();
}
}
void PreloadSource::clear()
{
{
QMutexLocker l(&_fileListMutex);
_filesToLoad.clear();
}
if (!isRunning()) {
start(QThread::LowPriority);
} else {
_condition.wakeOne();
}
//waitFinished.wait(&preloadMutex);
}
inline double PreloadSource::polyArea(const QVector<QPointF>& p) const
{
//https://en.wikipedia.org/wiki/Polygon#Area_and_centroid
const int n = p.size();
double area = 0.0;
for (int i=0; i<n; i++)
{
area += p[i].x()*p[(i+1)%n].y() - p[(i+1)%n].x()*p[i].y();
}
if (area < 0)
return -0.5*area;
else
return 0.5*area;
}
inline double PreloadSource::polyAreaCV(const std::vector<cv::Point> &p,
const int from,
const int to) const
{
//https://en.wikipedia.org/wiki/Polygon#Area_and_centroid
const int n = p.size();
const int localTo = to < 0 ? n : qMin(n, to);
double area = 0.0;
if (from < to) {
for (int i=from; i<localTo; i++)
{
const int toInd = i == localTo -1 ? from : i+1;
area += p[i].x*p[toInd].y - p[toInd].x*p[i].y;
}
} else {
// with walk over 0
for (int i=0; i<=from; i++)
{
const int toInd = i == from ? to : i+1;
area += p[i].x*p[toInd].y - p[toInd].x*p[i].y;
}
for (int i=localTo; i<n; i++)
{
const int toInd = i == n-1 ? 0 : i+1;
area += p[i].x*p[toInd].y - p[toInd].x*p[i].y;
}
}
if (area < 0)
return -0.5*area;
else
return 0.5*area;
}
void PreloadSource::getSumOfRectangleSampling(const QVector<cv::Point2f> pts,
const cv::Mat& image,
double& sumValues,
long int& numPixels)
{
// assume single channel image
sumValues = 0;
numPixels = 0;
IplImage tmp = image;
assert(image.type() == CV_8U);
for (int i=0; i<4; i++)
{
CvLineIterator iterator;
const int count = cvInitLineIterator( &tmp, pts[i], pts[(i+1)%4], &iterator, 4);
for( int j = 0; j < count; j++ ){
sumValues += iterator.ptr[0];
CV_NEXT_LINE_POINT(iterator);
}
/* for float images:
cv::Point2f diff = pts[(i+1)%4] - pts[i];
const int count= cv::norm(diff);
for (int j=0; j<count; j++) {
cv::Point2f p2(pts[i].x + diff.x * float(j) / (float)count,
pts[i].y + diff.y * float(j) / (float)count);
cv::Point2i p((int)p2.x, (int)p2.y);
if (p.x >=0 && p.y >=0
&& p.x < image.size().width
&& p.y < image.size().height)
{
sumValues += image.at<float>(p.y, p.x);
}
} */
numPixels += cv::norm(pts[i] - pts[(i+1)%4]);
}
}
void PreloadSource::getSumOfRectangleSamplingF(const QVector<cv::Point2f> pts,
const cv::Mat& image,
double& sumValues,
long int& numPixels)
{
// assume single channel image
sumValues = 0;
numPixels = 0;
assert(image.type() == CV_32F);
for (int i=0; i<4; i++)
{
// for float images:
cv::Point2f diff = pts[(i+1)%4] - pts[i];
const int count= cv::norm(diff);
for (int j=0; j<count; j++) {
cv::Point2f p2(pts[i].x + diff.x * float(j) / (float)count,
pts[i].y + diff.y * float(j) / (float)count);
cv::Point2i p((int)p2.x, (int)p2.y);
if (p.x >=0 && p.y >=0
&& p.x < image.size().width
&& p.y < image.size().height)
{
const float val = image.at<float>(p.y, p.x);
sumValues += val*0 == 0? val : 0;
}
}
numPixels += cv::norm(pts[i] - pts[(i+1)%4]);
}
}
struct QPairFirstComparer
{
template<typename T1, typename T2>
bool operator()(const QPair<T1,T2> & a, const QPair<T1,T2> & b) const
{
return a.first > b.first;
}
};
void PreloadSource::preSelectFast(const cv::Mat& image,
const std::vector<std::vector< cv::Point > >& contours,
const std::vector<cv::Vec4i>& hierarchy,
const QVector<cv::RotatedRect>& rects,
const QVector<QVector<cv::Point2f> >& points,
QVector<bool>& valid) const
{
const float imageArea = image.rows * image.cols;
int num =0;
#pragma omp parallel for
for (unsigned int i=0; i< contours.size(); i++)
{
if (!valid[i])
continue;
int level = 0;
int parent = hierarchy[i][3];
while (parent >= 0)
{
level++;
parent = hierarchy[parent][3];
if (level > _maxHierarchyLevel) {
break;
}
}
if (rects[i].size.area() / (double) imageArea < _minArea
/*|| ratio > _maxAspect do not check aspect since we may want to split this rectangle
which might result in a better aspect */
//a < _minAreaWithinImage
|| level > _maxHierarchyLevel)
{
#pragma omp critical
valid[i] = false;
num++;
}
}
}
#ifdef DEBUG
static int extractNum = 0;
#endif
inline float extractRectangleValueSum(const cv::Mat& image,
const cv::RotatedRect rect)
{
if (rect.size.area() <20) {
return 0;
}
cv::Mat M, rotated, cropped;
float angle = rect.angle;
cv::Size rect_size = rect.size;
// thanks to http://felix.abecassis.me/2011/10/opencv-rotation-deskewing/
if (rect.angle < -45.) {
angle += 90.0;
qSwap(rect_size.width, rect_size.height);
}
M = getRotationMatrix2D(rect.center, angle, 1.0);
// perform the affine transformation
cv::warpAffine(image, rotated, M, image.size(), cv::INTER_LINEAR);
// crop the resulting image
cv::getRectSubPix(rotated, rect_size, rect.center, cropped);
const float out = cv::sum(cropped)[0];
#ifdef DEBUG
cv::imwrite(QString("test_extract%1_e%2.png").arg(extractNum++).arg(out/(float)rect.size.area()).toLocal8Bit().toStdString(), cropped);
#endif
return out;
}
void PreloadSource::bestRectangles(const cv::Mat& imageOrig,
const cv::Mat& boundaries,
const cv::Vec3i &thresholds,
const std::vector<std::vector< cv::Point > >& contours,
const std::vector<cv::Vec4i>& hierarchy,
const QVector<cv::RotatedRect>& rects,
const QVector<QVector<cv::Point2f> >& points,
QVector<double>& energies,
QVector<bool> &out)
{
_minFrameEnergy = 0.75;
_maxAreaFractionEnergy = 0.5;
_e1 = 0.1;
_e2 = 5.0;
_e3 = 1;
_e4 = 10;
#ifdef DEBUG
cv::Mat bdCopy;
boundaries.copyTo(bdCopy);
cv::cvtColor(bdCopy, bdCopy, CV_GRAY2BGR);
#endif
const int imSize = 500;
cv::Mat smallImage;
// downsample very small, so it is feasible to extract whole rectangles
// there's usually a light peak from the scanner background
// hence compute energy better = higher difference from background
// this will be done later, however, he we prepare the corresponding image
float scale;
if (imageOrig.size().width < imageOrig.size().height)
{
scale = (float) imageOrig.size().height / (float) imSize;
cv::resize(imageOrig, smallImage, cv::Size(int(float(imageOrig.size().width)/scale),
imSize));
}
else
{
scale = (float) imageOrig.size().width / (float) imSize;
cv::resize(imageOrig, smallImage, cv::Size(imSize,
int(float(imageOrig.size().height)/scale)));
}
cv::Mat diffFromScannerBackground,
difference = cv::Mat::zeros(smallImage.size(), CV_8U) * 255.0,
differenceDark = cv::Mat::zeros(smallImage.size(), CV_8U) * 255.0,
tmp, backgroundEgdes;
std::vector<cv::Mat> channels(3);
cv::split(smallImage, channels);
for (int i=0; i<3; i++) {
diffFromScannerBackground = cv::Mat::ones(smallImage.size(), CV_8U) * thresholds[i];
// bright scanner background
cv::threshold(channels[i], tmp, thresholds[i], 0, cv::THRESH_TRUNC);
cv::absdiff(diffFromScannerBackground, tmp, diffFromScannerBackground);
difference = cv::max(difference, diffFromScannerBackground);
// black scanner background
diffFromScannerBackground = cv::Mat::ones(smallImage.size(), CV_8U) * 15; //thresholds[_histPositions][i];
cv::absdiff(diffFromScannerBackground, channels[i], diffFromScannerBackground);
differenceDark = cv::max(differenceDark, diffFromScannerBackground);
}
// * don't use the black values for now - many images seem dark and that would exclude quite some of them
// difference = cv::min(differenceDark, difference);
cv::medianBlur(difference, difference, 3);
const int threshTo = 255;
difference = difference * (threshTo / 50);
const int s = 3;
cv::Mat kernel = cv::getStructuringElement(0, cv::Size(s,s));
cv::Canny(difference, backgroundEgdes, 150, 300);
const int w = backgroundEgdes.size().width-1, h = backgroundEgdes.size().height-1;
cv::line(backgroundEgdes, cv::Point2i(0,0), cv::Point2i(w,0), cv::Scalar(255,2555,255,255));
cv::line(backgroundEgdes, cv::Point2i(0,0), cv::Point2i(0,h), cv::Scalar(255,2555,255,255));
cv::line(backgroundEgdes, cv::Point2i(w,0), cv::Point2i(w,h), cv::Scalar(255,2555,255,255));
cv::line(backgroundEgdes, cv::Point2i(0,h), cv::Point2i(w,h), cv::Scalar(255,2555,255,255));
cv::dilate(backgroundEgdes, backgroundEgdes, kernel);
cv::GaussianBlur(backgroundEgdes, backgroundEgdes, cv::Size(0,0), 3.5);
cv::threshold(difference, difference, 127, 255, cv::THRESH_BINARY);
#ifdef DEBUG_TEST
cv::imwrite(QString("test_smallCanny%1.png").arg(tmpVal2).toLocal8Bit().toStdString(), backgroundEgdes);
cv::imwrite(QString("test_smallImDiff%1.png").arg(tmpVal2++).toLocal8Bit().toStdString(), difference);
#endif
#ifdef DEBUG
cv::imwrite(QString("test_smallIm%1.png").arg(tmpVal).toLocal8Bit().toStdString(), smallImage);
#endif
int n = contours.size();
QVector<int> selected(n); // 0 not considered, 1 selected, -1 disregarded
for(int i=0; i< n; i++)
{
selected[i] = out[i] ? 0 : -1;
}
energies.resize(n);
#define INVALIDATE(i) \
out[i] = false; \
selected[i] = -1;
QVector<QPair<float, int> > energiesAndIdexes(n); // for sorting with backmap
const double imageArea = imageOrig.rows * imageOrig.cols;
// compute energy
#pragma omp parallel for
for (int i=0; i<n; i++)
{
if (!_stopped && !_abortAndLoadNext)
{
n = contours.size();
bool found = false;
double energy = 0;
if (out[i]) {
// ensure maximum aspect ratio ?
const double ratio = rects[i].size.width > rects[i].size.height ?
(double)rects[i].size.width / (double)rects[i].size.height :
(double)rects[i].size.height / (double)rects[i].size.width;
// rectangle is too much outside of the image?
QPolygonF pIm, pRect;
const float x = imageOrig.size().width, y = imageOrig.size().height;
pRect << CvPt2QtPt(points[i][0]) << CvPt2QtPt(points[i][1])
<< CvPt2QtPt(points[i][2]) << CvPt2QtPt(points[i][3]);
pIm << CvPt2QtPt(cv::Point2f(0,0)) << CvPt2QtPt(cv::Point2f(x,0))
<< CvPt2QtPt(cv::Point2f(x,y)) << CvPt2QtPt(cv::Point2f(0,y));
QPolygonF intersected = pRect.intersected(pIm);
const double withinImageArea = polyArea(intersected);
const double a = withinImageArea / rects[i].size.area();
int level = 0;
int parent = hierarchy[i][3];
while (parent >= 0)
{
level++;
parent = hierarchy[parent][3];
}
// all together + min Area
if (rects[i].size.area() / (double) imageArea > _minArea
&& rects[i].size.area() / (double) imageArea < _maxArea
&& ratio <= _maxAspect
&& out[i]
&& a > _minAreaWithinImage
&& level <= _maxHierarchyLevel)
{
const cv::RotatedRect& r = rects[i];
const cv::RotatedRect rectSmall(cv::Point2f(float(r.center.x)/scale,
float(r.center.y)/scale),
cv::Size2f(float(r.size.width)/scale,
float(r.size.height)/scale),
r.angle);
const double frameEnergy = extractRectangleValueSum(difference, rectSmall)
/ double((double)threshTo * (double)rectSmall.size.area());
if (frameEnergy > _minFrameEnergy)
{
double sumValues, backgroundSum;
long int numPixels, backgroundNum;
getSumOfRectangleSampling(points[i], boundaries, sumValues, numPixels);
getSumOfRectangleSampling(points[i], backgroundEgdes, backgroundSum, backgroundNum);
#ifdef DEBUG
QVector<cv::Point2f> pts(4);
rectSmall.points(pts.data());
for (int k=0; k<4; k++) {
cv::line(bdCopy, points[i][k], points[i][(k+1)%4], cv::Scalar(0,0,255));
cv::line(smallImage, pts[k], pts[(k+1)%4], cv::Scalar(0,0,255));
}
#endif
// best area is around half of the image
// energy: 1. sampling values along the field
// 2. area (highest value is half the image)
// 3. aspect - 7:5 (between 3:2 and 4:3)
// 4. monotony inside rect (used above the reject background rectangles)
const double areaEnergy = (imageArea*_maxAreaFractionEnergy - qAbs(double(imageArea*_maxAreaFractionEnergy - rects[i].size.area()))) / (imageArea*_maxAreaFractionEnergy);
const double borderEnergy = sumValues / (double)numPixels / 255.0;
const double borderEnergyBackground = backgroundSum / (double)backgroundNum / 255.0;
const double aspectEnergy = (7.0/5.0 - qAbs(7.0/5.0 - ratio)) * 5.0 / 7.0;
energy = _e1 * areaEnergy * areaEnergy
+ _e2 * qPow(borderEnergy, 2.0)
+ _e2 * qPow(borderEnergyBackground, 2.0)
+ _e3 * aspectEnergy * aspectEnergy
+ _e4 * frameEnergy * frameEnergy;
/*
double tmp = areaEnergy, tmp2 = aspectEnergy, pow =2.0;
energy = _e1 * qPow(0.2*qRound(tmp*5.0), pow)
+ _e2 * qPow(qMin(borderEnergy, frameEnergy ), pow)
+ _e4 * qPow(0.2*qRound(tmp2*5.0), pow);*/
found = true;
}
}
}
#pragma omp critical
if (found) {
energiesAndIdexes[i] = QPair<float,int>(energy, i);
energies[i] = energy;
}
else
{
energiesAndIdexes[i] = QPair<float,int>(0, i);
INVALIDATE(i);
}
}
}
// Ordering ascending
qSort(energiesAndIdexes.begin(), energiesAndIdexes.end(), QPairFirstComparer());
// now choose: best first, then disable all parents and kids. search next
int choice = 0;
while (choice >= 0 && !_stopped && !_abortAndLoadNext)
{
choice = -1;
// search next best rectangle
for (int i=0; i<n; i++)
{
if (selected[energiesAndIdexes[i].second] == 0) {
choice = energiesAndIdexes[i].second;
break;
}
}
if (choice <=0) {
break;
}
selected[choice] = 1;
out[choice] = true;
// hierarchy[i]: next(0), prev(1), firstKid(2), parent(3)
// look for parents
int parent = hierarchy[choice][3];
while(parent >=0)
{
INVALIDATE(parent);
parent = hierarchy[parent][3];
}
// now stack kids and work through
QList<int> kids;
if (hierarchy[choice][2]>=0)
{
kids.append(hierarchy[choice][2]);
}
while (kids.length()>0)
{
const int current = kids.first();
kids.pop_front(); // remove that element
if (current < 0)// || selected[current] <0)
{
continue;
}
INVALIDATE(current);
// append first of that current kid
kids.append(hierarchy[current][2]);
// append next neighbour of that current kid
kids.append(hierarchy[current][0]);
}
}
#ifdef DEBUG
cv::imwrite(QString("test_rectsSampling%1.png").arg(tmpVal).toLocal8Bit().toStdString(), bdCopy);
cv::imwrite(QString("test_rectsSamplingSmall%1.png").arg(tmpVal).toLocal8Bit().toStdString(), smallImage);
#endif
}
inline QVector<float> computeContourLength(const std::vector<cv::Point>& contour)
{
const int m = contour.size();
// compute sum of contour lengths
QVector<float> length(m+1);
for (int j=1; j<m; j++)
{
length[j] = length[j-1] + cv::norm(contour[j] - contour[j-1]);
}
length[m] = length[m-1] + cv::norm(contour[0] - contour[m-1]);
return length;
}
inline void addSubRectangles(const int i,
const int i1,
const int i2,
std::vector<std::vector< cv::Point > >& contours,
QVector<cv::RotatedRect>& rects,
std::vector<cv::Vec4i>& hierarchy,
QVector<QVector<cv::Point2f> >& points,
QVector<bool>& valid)
{
// ok here it seems there is a possibility such that two rectangles are merged into one
// now add these two sub-rectangles/contours and make sure hierarchy is just correct
// (however the "inside" relation won't be correct)
// hierarchy[i]: next(0), prev(1), firstKid(2), parent(3)
const int l= contours.size();
std::vector< cv::Point > c1, c2;
c1.insert(c1.end(), contours[i].begin(), contours[i].begin() + (i1==(int)contours[i].size()-1? i1: i1+1));
c1.insert(c1.end(), contours[i].begin() + i2, contours[i].end()); //(i2==0? i2 : i2-1), contours[i].end() );
c2.insert(c2.end(),
contours[i].begin() + i1, //(i1==0?i1:i1-1),
contours[i].begin() + (i2==(int)contours[i].size()-1? i2: i2+1));
contours.push_back(c1);
contours.push_back(c2);
rects.push_back(minAreaRect (c1));
rects.push_back(minAreaRect (c2));
valid.push_back(true);
valid.push_back(true);
points.resize(l+2);
points[l].resize(4);
points[l+1].resize(4);
rects[l].points(points[l].data());
rects[l+1].points(points[+1l].data());
hierarchy.push_back(cv::Vec4i(l+1, -1, -1, i));
if (hierarchy[i][2] >= 0)
{
const int kid = hierarchy[i][2];
hierarchy[kid][1] = l+1;
hierarchy.push_back(cv::Vec4i(kid, l, -1, i));
}
else
{
hierarchy.push_back(cv::Vec4i(-1, l, -1, i));
}
hierarchy[i][2] = l;
}
void PreloadSource::rectangleOverlap(std::vector<std::vector< cv::Point > >& contours,
QVector<cv::RotatedRect>& rects,
std::vector<cv::Vec4i>& hierarchy,
QVector<QVector<cv::Point2f> >& points,
QVector<QVector<float> > &sumOfContourLengths,
QVector<bool>& valid) const
{
// find rectangles which (slightly) overlap at two respective corners
// also look for single large convexity defects
const int n = contours.size();
for (int i=0; i<n && !_stopped && !_abortAndLoadNext; i++)
{
if (!valid[i])
continue;
std::vector<cv::Vec4i> defects;
std::vector<int> hull;
cv::convexHull(contours[i], hull);
if (hull.size() < 4)
{
continue;
}
cv::convexityDefects(contours[i], hull, defects);
if (defects.size() > 1)
{
// find best defect pair, which means closest distance to each other
// and long enough contour length in between
const double rectLength = 2*rects[i].size.width + 2*rects[i].size.height;
const int m = contours[i].size();
// now compare defects, but only if contour in between is long enough
for (unsigned int j=0; j < defects.size(); j++)
{
if (sumOfContourLengths[i].size() < (int)contours[i].size()) {
sumOfContourLengths[i] = computeContourLength(contours[i]);
}
const QVector<float>& length = sumOfContourLengths[i];
// getMaxDefect position
const double maxDist = qMin(rects[i].size.width,
rects[i].size.height)/2.0;
const int bestN = 5;
QVector<double> maxDistFound (bestN, 0);
QVector<QPair<int, int> > bestPairs(bestN);
/**
* compare to closest point on rectangle. if it is only (say half)
* the length away (-> to the opposite side) then add splitted
*/
//const double minLengthInBetween = 40.0;
//const double maxAreaPerLength = 5.0;
const int ind = defects[j][2];
for (unsigned int j=0; j<contours[i].size(); j++)
{
const int i1 = qMin(ind, (int)j), i2 = qMax(ind, (int)j);
const double d1 = length[i2] - length[i1];
const double d2 = length[i1] + (length[m] - length[i1]);
double shortestDist = qMin(d1, d2);
if (cv::norm(contours[i][i1] - contours[i][i2]) < maxDist
&& shortestDist > rectLength * _splitMinLengthFrac)
{
// energy consists of length fraction and distance to rectangle fraction
// so defects with rather equal lengths towards both circles
// and a high defect values are preferred
const double distL1P1 = distancePointLine(points[i][0], points[i][1], contours[i][i1]);
const double distL1P2 = distancePointLine(points[i][2], points[i][3], contours[i][i1]);
const double distL2P1 = distancePointLine(points[i][1], points[i][2], contours[i][i1]);
const double distL2P2 = distancePointLine(points[i][0], points[i][3], contours[i][i1]);
double minRectDist = qMin(distL1P1, qMin(distL1P2, qMin(distL2P1, distL2P2)));
// try to distribute around the rectangle
const double energy = (shortestDist - length.last() / (double)bestN)
+ 4 * minRectDist * rectLength;
int choice = bestN;
for (int k=bestN-1; k>=0; k--) {
if (energy < maxDistFound[k] )
{
break;
}
choice = k;
}
if (choice < bestN) {
bestPairs[choice] = QPair<int, int>(i1, i2);
maxDistFound[choice] = energy;
}
}
}
for (int k=0; k<bestN; k++) {
if (maxDistFound[k] > 0) {
addSubRectangles(i, bestPairs[k].first, bestPairs[k].second, contours, rects, hierarchy, points, valid);
}
}
/**
* the following seems to work ok, however it is somewhat limited
* since it is only tried to compare to neighboring defects
*/
// compare to other defects
for (unsigned int k=j+1; k < defects.size(); k++)
{
const int i1 = qMin(defects[j][2], defects[k][2]);
const int i2 = qMax(defects[j][2], defects[k][2]);
double shortestDist = qMin(length[i2] - length[i1],
length[i1] + (length[m] - length[i1]));
if (shortestDist > rectLength * _splitMinLengthFrac
&& cv::norm(contours[i][i1] - contours[i][i2]) < maxDist)
{
// now save this
// found best defect pair
// now test how close that pair is and if it is rather on the diagonal of the bounding rect
const double diagLength = qSqrt(rects[i].size.width * rects[i].size.width
+ rects[i].size.height * rects[i].size.height);
const double distD1P1 = distancePointLine(points[i][0], points[i][2], contours[i][i1]);
const double distD1P2 = distancePointLine(points[i][0], points[i][2], contours[i][i2]);
const double distD2P1 = distancePointLine(points[i][1], points[i][3], contours[i][i1]);
const double distD2P2 = distancePointLine(points[i][1], points[i][3], contours[i][i2]);
const double sumDiagDist = qMin( distD1P1 + distD1P2, distD2P1 + distD2P2);
if (sumDiagDist < _splitMaxOffsetFrac * diagLength)
{
// also make sure this is far enough away from corner
double distanceCorner = 1e10;
for (int l=0; l<4; l++)
{
distanceCorner = qMin(distanceCorner, cv::norm(points[i][l]));
}
if (distanceCorner > _splitMinCornerDist * diagLength) {
addSubRectangles(i, i1, i2, contours, rects, hierarchy, points, valid);
}
}
}
}
}
}
}
}
/*
void PreloadSource::computeContourLengths(const std::vector<std::vector< cv::Point > >& contours,
const QVector<bool>& valid,
QVector<QVector<float> >& sumOfContourLengths)
{
const int n = contours.size();
sumOfContourLengths.resize(n);
#pragma omp parallel for
for (int i=0; i<n; i++)
{
if (valid[i])
{
const int m = contours[i].size();
// compute sum of contour lengths
QVector<float> length(m+1);
for (int j=1; j<m; j++)
{
length[j] = length[j-1] + cv::norm(contours[i][j] - contours[i][j-1]);
}
length[m] = length[m-1] + cv::norm(contours[i][0] - contours[i][m-1]);
#pragma omp critical
sumOfContourLengths[i] = length;
}
}
}
*/
void PreloadSource::removeFilaments(QVector<QVector<float> >& sumOfContourLengths,
const QVector<bool> & valid,
std::vector<std::vector< cv::Point > >& contours)
{
/*
*I did remove this since there's quite some danger correct images might get removed, as well
*/
const int n = contours.size();
#pragma omp parallel for
for (int i=0; i<n; i++)
{
if (!valid[i]) {
continue;
}
if (sumOfContourLengths[i].size() < (int)contours[i].size()) {
sumOfContourLengths[i] = computeContourLength(contours[i]);
}
const int m = contours[i].size();
QVector<bool> valid(m, true);
// compute all distances
for (int j=1; j<m; j++)
{
if (valid[j])
{
for (int k=j+1; k<m; k++)
{
if (valid[k])
{
double dist = cv::norm(contours[i][j] - contours[i][k]);
const double maxDist = 10.0;
const double minLengthInBetween = 40.0;
const double maxAreaPerLength = 1.0;
if (dist <= maxDist) {
const double cLength1 = sumOfContourLengths[i][k] -
sumOfContourLengths[i][j];
if (cLength1 >= minLengthInBetween) {
// now compute area of this part polygon.
// if it is very small, then add this pair to removal.
// sample both directions
const double areaPerLength1 =
polyAreaCV(contours[i], j, k ) / cLength1;
if (areaPerLength1 < maxAreaPerLength)
{
// prevent computing this, again
for (int v=j+1; v<k; v++) valid[v] = false;
}
}
const double cLength2 = sumOfContourLengths[i][j] +
(sumOfContourLengths[i].last() - sumOfContourLengths[i][j]);
if (cLength2 >= minLengthInBetween)
{
const double areaPerLength2 =
polyAreaCV(contours[i], k, j ) / cLength2;
if (areaPerLength2 < maxAreaPerLength)
{
// prevent computing this, again
for (int v=k+1; v<m; v++) valid[v] = false;
for (int v=0; v<j; v++) valid[v] = false;
}
}
}
}
}
}
}
#pragma omp critical
{
for (int j=m-1; j>=0; j--) {
if (!valid[j]) {
contours[i].erase(contours[i].begin() + j);
}
}
}
}
/**/
}
double PreloadSource::distancePointLine(const cv::Point2f& LP1,
const cv::Point2f& LP2,
const cv::Point2f& p) const
{
const double dY = LP2.y - LP1.y;
const double dX = LP2.x - LP1.y;
const double ret = qAbs ( dY * p.x - dX * p.y + LP2.x * LP1.y - LP2.y * LP1.x) / qSqrt(dY * dY + dX * dX);
return ret;
}
cv::Mat PreloadSource::loadAndShrink(const QString& filename)
{
cv::Mat image = cv::imread(filename.toLocal8Bit().data(),
CV_LOAD_IMAGE_ANYDEPTH | CV_LOAD_IMAGE_ANYCOLOR);
if (image.channels() == 1) {
cv::cvtColor(image, image, CV_GRAY2BGR);
}
if (image.channels() == 4) {
cv::cvtColor(image, image, CV_BGRA2BGR);
}
if (image.depth() == CV_16S
|| image.depth() == CV_16U)
{
qDebug() << "found 16bit image";
cv::normalize(image,image,0,255,cv::NORM_MINMAX);
image.convertTo(image, CV_8U);
}
return image;
}
void PreloadSource::newTargets(SourceFilePtr &source, cv::Mat& alreadyLoaded, bool& locked)
{
source->targets.clear();
_imSize = 2000;
// prevent two big chunks of memory at the same time -> see copytargets
if (alreadyLoaded.empty())
{
alreadyLoaded = loadAndShrink(source->source.absoluteFilePath());
}
cv::Mat image = alreadyLoaded;
if (image.size().width == 0
|| image.size().height == 0
|| image.empty())
{
source->error = true;
return;
}
float scale;
cv::Mat im;
if (image.size().width < image.size().height)
{
scale = (float) image.size().height / (float) _imSize;
cv::resize(image, im, cv::Size(int(float(image.size().width)/scale),
_imSize));
}
else
{
scale = (float) image.size().width / (float) _imSize;
cv::resize(image, im, cv::Size(_imSize,
int(float(image.size().height)/scale)));
}
CopyTargets::memorySaveMutex.unlock();
locked = false;
ABORT_COMPUTATION();
int bestThresh;
std::vector<cv::Vec3i> thresholds = getThresholds(image, bestThresh);
std::vector<cv::RotatedRect> rectDetection = extractRectangles(thresholds, image);
cv::cvtColor(im, image, cv::COLOR_BGR2GRAY);
cv::medianBlur(image, image, 5);
image.convertTo(image, CV_32F);
if (image.cols == 0) {
return;
}
cv::Mat overallMask, out;
const int levels = _levels;
overallMask = cv::Mat::zeros(image.rows, image.cols, image.type());
for (int i=0; i < levels; i++)
{
ABORT_COMPUTATION()
cv::pyrDown(image, image);
cv::Mat mu, mu2, sigma;
cv::blur(image, mu, cv::Size(3,3));
cv::Mat tmp(image.rows, image.cols, image.type());
cv::multiply(image, image, tmp);
cv::blur(tmp, mu2, cv::Size(3,3));
cv::multiply(mu, mu, tmp);
cv::sqrt(mu2 - tmp, sigma);
cv::resize(sigma, tmp, cv::Size(overallMask.cols, overallMask.rows));
overallMask += tmp;
}
image.release();
// remove NaN and INF
for(int row = 0; row < overallMask.rows; ++row) {
float* p = overallMask.ptr<float>(row);
for(int col = 0; col < overallMask.cols; col++) {
if (p[col] * 0.0f != 0.0f) {
p[col] = 0.0f;
}
}
}
ABORT_COMPUTATION()
const int s = 5;
cv::Mat kernel = cv::getStructuringElement(0, cv::Size(s,s));
cv::morphologyEx(overallMask, overallMask, cv::MORPH_OPEN, kernel);
cv::erode(overallMask, overallMask, kernel);
cv::threshold(overallMask, out, _threshold*levels, 1, CV_THRESH_BINARY);
out.convertTo(out, CV_8U);
ABORT_COMPUTATION()
#ifdef DEBUG
// save mask
QString fn = source->source.baseName() + "_mask.png";
cv::imwrite(fn.toLocal8Bit().toStdString(), overallMask);
#endif
std::vector<std::vector< cv::Point > > contours;
std::vector<cv::Vec4i> hierarchy;
cv::findContours( out, contours, hierarchy,
CV_RETR_TREE , //CV_RETR_LIST,
CV_CHAIN_APPROX_TC89_KCOS //CV_CHAIN_APPROX_NONE
);
ABORT_COMPUTATION()
if (rectDetection.size() > 0)
{
const int start = contours.size() ;
contours.resize(start + rectDetection.size());
hierarchy.resize(start + rectDetection.size());
cv::Vec4i hi;
hi[0] = hi[1] = hi[2] = hi[3] = -1;
for (unsigned int i=0; i< rectDetection.size(); i++)
{
ABORT_COMPUTATION()
hierarchy[i+start] = hi;
contours[i+start].resize(4);
std::vector<cv::Point2f> pts(4);
rectDetection[i].size.width /= scale;
rectDetection[i].size.height /= scale;
rectDetection[i].center.x /= scale;
rectDetection[i].center.y /= scale;
rectDetection[i].points(pts.data());
for (int j=0; j<4; j++) {
contours[i+start][j] = pts[j];
}
}
}
// compute needed data
QVector<cv::RotatedRect> rects(contours.size());
QVector<QVector<cv::Point2f> > points(contours.size());
QVector<bool> keep(contours.size());
keep.fill(true);
QVector<double> energies;
for (int j=0; j< (int)contours.size(); j++) {
rects[j] = minAreaRect (contours[j]);
points[j].resize(4);
rects[j].points(points[j].data());
}
ABORT_COMPUTATION()
PreloadSource::preSelectFast(overallMask, contours, hierarchy, rects, points, keep);
ABORT_COMPUTATION()
#ifdef DEBUG
/// Draw contours
///
cv::RNG rng(12345);
cv::Mat drawing = cv::Mat::zeros( overallMask.size(), CV_8UC3 );
for( unsigned int i = 0; i< contours.size(); i++ )
{
if (hierarchy[i][3] <0) {
cv::Scalar color = cv::Scalar( rng.uniform(0, 255), rng.uniform(0,255), rng.uniform(0,255) );
cv::drawContours( drawing, contours, i, color, 2, 8, hierarchy, _maxHierarchyLevel);
}
}
fn = source->source.baseName() + "_conts.png";
cv::imwrite(fn.toLocal8Bit().toStdString(), drawing);
#endif
QVector<QVector<float> > sumOfContourLengths(contours.size());
//computeContourLengths(contours, sumOfContourLengths);
//removeFilaments(sumOfContourLengths, keep, contours);
// detect overlapping rectangles and split others, as well
rectangleOverlap(contours, rects, hierarchy, points, sumOfContourLengths, keep);
ABORT_COMPUTATION()
// select best rectangle based on hierachy structure
//cv::GaussianBlur(overallMask, overallMask, cv::Size(0,0), 3.5);
// do it again, find contours seems to destroy out
// cv::threshold(overallMask, out, _threshold*levels, 255, CV_THRESH_BINARY);
out = overallMask * int(_threshold * float(levels));
cv::GaussianBlur(out, out, cv::Size(0,0), 5.0);
out.convertTo(out, CV_8U);
#ifdef DEBUG
fn = source->source.baseName() + "_blur.png";
cv::imwrite(fn.toLocal8Bit().toStdString(), out);
#endif
ABORT_COMPUTATION()
bestRectangles(im,
out,
thresholds[bestThresh>=0? bestThresh : 0],
contours,
hierarchy,
rects,
points,
energies,
keep);
ABORT_COMPUTATION()
#ifdef DEBUG
/// Draw contours
///
drawing = cv::Mat::zeros( overallMask.size(), CV_8UC3 );
for( unsigned int i = 0; i< contours.size(); i++ )
{
if (hierarchy[i][3] <0) {
cv::Scalar color = cv::Scalar( rng.uniform(0, 255), rng.uniform(0,255), rng.uniform(0,255) );
cv::drawContours( drawing, contours, i, color, 2, 8, hierarchy, _maxHierarchyLevel);
}
}
fn = source->source.baseName() + "_conts_noFilaments.png";
cv::imwrite(fn.toLocal8Bit().toStdString(), drawing);
#endif
// now extract rectangles and remove remaining overlaps
// sort by energy. this means a rectangle cannot lose and the winner loses itself.
QVector<QPair<double, int> > backmap(contours.size());
for (unsigned int i=0; i<contours.size(); i++)
{
backmap[i] = QPair<double, int>(energies[i], i);
}
qSort(backmap.begin(), backmap.end(), QPairFirstComparer());
#ifdef DEBUG_TEST
cv::Mat visualize;
im.copyTo(visualize);
#endif
cv::GaussianBlur(overallMask, overallMask, cv::Size(0,0), 5.0);
for (int j=0; j< (int)contours.size(); j++) {
ABORT_COMPUTATION()
const int ind1 = backmap[j].second;
cv::RotatedRect& rect = rects[ind1];
if (keep[ind1]) {
// find overlapping rectangles - disregard the one with less energy
QPolygonF p;
p << CvPt2QtPt(points[ind1][0]) << CvPt2QtPt(points[ind1][1])
<< CvPt2QtPt(points[ind1][2]) << CvPt2QtPt(points[ind1][3]);
for (int i=j+1; i<(int)contours.size(); i++)
{
const int ind2 = backmap[i].second;
if (keep[ind2])
{
QPolygonF p2;
p2 << CvPt2QtPt(points[ind2][0]) << CvPt2QtPt(points[ind2][1])
<< CvPt2QtPt(points[ind2][2]) << CvPt2QtPt(points[ind2][3]);
const QPolygonF intersected = p.intersected(p2);
const double area = polyArea(intersected);
const double a1 = area / (rects[ind1].size.area());
const double a2 = area / (rects[ind2].size.area());
if (a1 > _maxOverlap || a2 > _maxOverlap)
{
keep[ind2] = false;
}
}
}
#ifdef DEBUG_TEST
rect.points(points[ind1].data());
for (int i=0; i<4; i++) {
cv::line(visualize, points[ind1][i], points[ind1][(i+1)%4], cv::Scalar(0,0,255), 3);
}
#endif
// optimize rectangle
optimizeRectangle(overallMask, rect);
rect.points(points[ind1].data());
#ifdef DEBUG_TEST
for (int i=0; i<4; i++) {
cv::line(visualize, points[ind1][i], points[ind1][(i+1)%4], cv::Scalar(0,255,255), 2);
}
#endif
const float scaleNew = scale / source->scale;
QVector<QPointF> corners(4);
corners[0] = QPointF(scaleNew*rect.size.width/2.0, scaleNew*rect.size.height/2.0);
corners[1] = QPointF(scaleNew*rect.size.width/2.0, -scaleNew*rect.size.height/2.0);
corners[2] = QPointF(-scaleNew*rect.size.width/2.0, -scaleNew*rect.size.height/2.0);
corners[3] = QPointF(-scaleNew*rect.size.width/2.0, scaleNew*rect.size.height/2.0);
TargetImagePtr target(new TargetImage);
target->boundary->setCorners(corners.data());
target->boundary->setTransform(target->boundary->transform().translate(scaleNew*rect.center.x, scaleNew*rect.center.y).rotate(rect.angle));
source->targets.push_back(target);
}
}
#ifdef DEBUG_TEST
QString fn2 = QString("%1_%2_%3_%4_" + source->source.baseName() + "_prodBorder_results.png").arg(_e1).arg(_e2).arg(_e3).arg(_e4);
cv::imwrite(fn2.toLocal8Bit().toStdString(), visualize);
#endif
}
void PreloadSource::updateMinMax(const int mi, const int ma)
{
_min = mi;
_max = ma;
}
void PreloadSource::stop()
{
_stopped = true;
}
QString PreloadSource::isCurrentlyLoading() const
{
return _currentFilename;
}
inline bool PreloadSource::testAbort(SourceFilePtr sourceFile)
{
if (_stopped || _abortAndLoadNext) {
sourceFile->imageOrig = QImage();
sourceFile->image.reset();
preloadMutex.unlock();
waitFinished.wakeAll();
_currentFilename = "";
return true;
}
return false;
}
#define TEST_ABORT(s) \
if (testAbort(s)) { \
if (locked) CopyTargets::memorySaveMutex.unlock(); \
locked = false; \
continue; \
}
void PreloadSource::run()
{
while (!_stopped) {
_abortAndLoadNext = false;
SourceFilePtr sourceFile;
int id;
{
QMutexLocker l(&_fileListMutex);
if (_filesToLoad.size() == 0) {
waitFinished.wakeOne();
_condition.wait(&_fileListMutex);
// continue: in the meantime it might have been stopped
// or set empty
continue;
}
sourceFile = _filesToLoad.first().first;
_current = sourceFile;
id = _filesToLoad.first().second;
_filesToLoad.pop_front();
}
preloadMutex.lock();
_currentFilename = sourceFile->source.absoluteFilePath();
qDebug() << "starting to load: " << _currentFilename;
QImage tmpIm = QImage(sourceFile->source.absoluteFilePath());
CopyTargets::memorySaveMutex.lock();
bool locked = true;
sourceFile->error = false;
sourceFile->imageOrig = tmpIm.convertToFormat(_qImageFmt);
TEST_ABORT(sourceFile);
cv::Mat loadedData;
if (!sourceFile->imageOrig.isNull() || sourceFile->imageOrig.width() == 0)
{
// try to load a 16 bit image with opencv
loadedData = loadAndShrink(sourceFile->source.absoluteFilePath());
if (loadedData.empty())
{
sourceFile->error = true;
}
else
{
// QImage has different rgb ordering
cv::Mat tmp;
cv::cvtColor(loadedData, tmp, CV_BGR2RGB);
sourceFile->imageOrig = QImage(tmp.data,
tmp.size().width,
tmp.size().height,
tmp.step[0],
QImage::Format_RGB888).convertToFormat(_qImageFmt);
}
}
TEST_ABORT(sourceFile);
if (!sourceFile->error)
{
// scale down original image
const int maxExtent = 2000;
sourceFile->scale = 1.0;
if (sourceFile->imageOrig.width() > maxExtent
&& sourceFile->imageOrig.width() >= sourceFile->imageOrig.height())
{
sourceFile->scale = (double) sourceFile->imageOrig.width() / (double)maxExtent;
sourceFile->imageOrig = sourceFile->imageOrig.scaled(QSize(maxExtent, maxExtent), Qt::KeepAspectRatio);
}
else if (sourceFile->imageOrig.height() > maxExtent
&& sourceFile->imageOrig.width() < sourceFile->imageOrig.height())
{
sourceFile->scale = (double) sourceFile->imageOrig.height() / (double)maxExtent;
sourceFile->imageOrig = sourceFile->imageOrig.scaled(QSize(maxExtent, maxExtent), Qt::KeepAspectRatio);
}
TEST_ABORT(sourceFile);
if (id >= 0)
{
newTargets(sourceFile, loadedData, locked); // if id < 0 then it means only load images
}
}
if (locked) CopyTargets::memorySaveMutex.unlock();
locked = false;
TEST_ABORT(sourceFile);
// remove from list - since we will have it commited in a couple of lines
#ifndef TEST_DEBUG
if (!_stopped)
{
emit doneLoading(sourceFile, id < 0 ? -id -1 : id);
}
#else
sourceFile->image = QGraphicsPixmapItemPtr();
sourceFile->imageOrig = QImage();
for (int i=0; i<sourceFile.targets.size(); i++) {
sourceFile->targets[i]->image = QImage();
}
sourceFile->targets.clear();
#endif
_currentFilename = "";
preloadMutex.unlock();
waitFinished.wakeAll();
}
}
/**
* advanced extracting here
*/
inline cv::Point2f lineIntersection(const cv::Vec2f& l1,
const cv::Vec2f& l2)
{
const float& theta1 = l1[1];
const float& rho1 = l1[0];
const float& theta2 = l2[1];
const float& rho2 = l2[0];
const float c1 = qCos(theta1);
const float c2 = qCos(theta2);
const float s1 = qSin(theta1);
const float s2 = qSin(theta2);
/* cross product l1 x l2 = intersection point (homogeneous)
[r1*s2 - r2*s1],
[c1*r2 - c2*r1],
[c1*s2 - c2*s1]] */
cv::Point2f out;
const float third = c1*s2 - c2*s1;
out.x = (rho1*s2 - rho2*s1)/third;
out.y = (c1*rho2 - c2*rho1)/third;
return out;
}
inline std::vector<cv::Point2f> rectanglePoints(const cv::Vec2f& l1,
const cv::Vec2f& l2,
const cv::Vec2f& l3,
const cv::Vec2f& l4)
{
std::vector<cv::Point2f> out(4);
out[0] = lineIntersection(l1, l2);
out[1] = lineIntersection(l2, l3);
out[2] = lineIntersection(l3, l4);
out[3] = lineIntersection(l4, l1);
return out;
}
std::vector<float> lineValues(const cv::Mat& image,
const cv::Point2f p1,
const cv::Point2f p2)
{
assert(image.type() == CV_32F);
cv::Point2f diff = p2 - p1;
const int count= cv::norm(diff);
std::vector<float> out(count);
for (int i=0; i<count; i++) {
cv::Point2f p(p1.x + diff.x * float(i) / (float)count,
p1.y + diff.y * float(i) / (float)count);
out[i] = image.at<float>((int)p.y, (int)p.x);
}
return out;
/*
const IplImage tmp = image;
CvLineIterator iterator;
const int count = cvInitLineIterator( &tmp, p1, p2, &iterator, connectivity);
std::vector<float> out(count);
for( int i = 0; i < count; i++ ){
out[i] = iterator.ptr[0];
CV_NEXT_LINE_POINT(iterator);
}
return out;*/
}
inline float dot(const cv::Vec2f& v1, const cv::Vec2f v2)
{
return v1[0]*v2[0] + v1[1]*v2[1];
}
std::vector<cv::RotatedRect> PreloadSource::extractRectangles(const std::vector<cv::Vec3i> &thresholds,
const cv::Mat& image)
{
_rectAngularThresh = M_PI * 2.0/180.0;
_rectMinLineDistFrac = 0.3;
_maxArea = 0.75;
const float minLineDist = qMin(image.size().width, image.size().height) * _rectMinLineDistFrac;
std::vector<cv::Vec2f> lines = extractLines(thresholds, image);
if (lines.size() >= 150)
{
const int n= lines.size();
std::vector<cv::Vec2f> lines2;
lines2.reserve(n/150);
for (unsigned int i=0; i<lines.size(); i+= (int) n /150)
{
lines2.push_back(lines[i]);
}
lines = lines2;
}
_rectMinDist = image.size().width > image.size().height ?
image.size().width/2000
: image.size().height/2000 ;
_rectMinDist = qMax(_rectMinDist, 2.0f);
// add the boundary lines, since images at the boundary
// will usually have no lines detected at that boundary
lines.reserve(lines.size() + 4);
lines.push_back(cv::Vec2f(0, 0.5 * M_PI));
lines.push_back(cv::Vec2f(0, 2.0));
lines.push_back(cv::Vec2f(image.size().width, 0.0));
lines.push_back(cv::Vec2f(image.size().height, 0.5 * M_PI));
const int n = lines.size();
QList<QPair<int, int> > parallelLines;
QList<float> averageAngle;
#pragma omp parallel for
for (int i=0; i<n; i++) {
const float& rho1 = lines[i][0];
const float& theta1 = lines[i][1];
for (int j= i+1; j < n; j++) {
const float& rho2 = lines[j][0];
const float& theta2 = lines[j][1];
const float angularDistance = qAbs(theta1 - theta2);
const float lineDist = qAbs(rho1 - rho2);
if (angularDistance < _rectAngularThresh
&& lineDist > minLineDist)
{
#pragma omp critical
{
parallelLines.append(QPair<int, int>(i,j));
averageAngle.append(0.5 * (theta1 + theta2));
}
}
}
}
/*
* look for pairs of parallel lines which meet the rectangle criteria
*/
const float imageArea = image.size().width * image.size().height;
std::vector<cv::RotatedRect> out;
std::vector<std::vector<cv::Point2f> > points;
QList<QPair<int, int> >::const_iterator it1(parallelLines.begin());
QList<float>::const_iterator avgAngle1(averageAngle.begin());
while (it1 != parallelLines.end() && !_stopped && !_abortAndLoadNext)
{
const QPair<int, int> parallelLines1 = *it1; it1++;
const float angle1 = *avgAngle1; avgAngle1++;
QList<QPair<int, int> >::const_iterator it2(it1);
QList<float>::const_iterator avgAngle2 (avgAngle1);
while(it2 != parallelLines.end())
{
const QPair<int, int> parallelLines2 = *it2; it2++;
const float angle2 = *avgAngle2; avgAngle2++;
// 90 degree
const float angularDist = qAbs(qAbs(angle1 - angle2) - 0.5 * M_PI);
if (angularDist < _rectAngularThresh)
{
// it is close to a rectangle
// estimate enclosing rectangle points
std::vector<cv::Point2f> pts = rectanglePoints(
lines[parallelLines1.first],
lines[parallelLines2.first],
lines[parallelLines1.second],
lines[parallelLines2.second]);
const cv::RotatedRect rect = cv::minAreaRect(pts);
const float aspect = rect.size.width > rect.size.height ?
rect.size.width / rect.size.height
: rect.size.height / rect.size.width;
assert(aspect >= 1.0);
const float relArea = rect.size.area() / imageArea;
if (aspect > _maxAspect || relArea < _minArea || relArea > _maxArea)
{
continue;
}
out.push_back(rect);
points.push_back(pts);
}
}
}
#ifdef DEBUG
cv::Mat im;
image.copyTo(im);
for (unsigned int i=0; i<lines.size(); i++)
{
float a = qCos(lines[i][1]);
float b = qSin(lines[i][1]);
float x0 = a*lines[i][0];
float y0 = b*lines[i][0];
int x1 = int(x0 + 10000*(-b));
int y1 = int(y0 + 10000*(a));
int x2 = int(x0 - 10000*(-b));
int y2 = int(y0 - 10000*(a));
cv::line( im, cv::Point(x1,y1), cv::Point(x2,y2), cv::Scalar(255,0,0), 2, 8 );
}
for (unsigned int i=0; i<out.size(); i++)
{
std::vector<cv::Point2f> pts(4);
out[i].points(pts.data());
cv::Scalar colorScalar = cv::Scalar( (i*7)%256, 128+ (i*19)%256, 220 + (i*11 )%256);
for( int j = 0; j < 4; j++ )
cv::line( im, pts[j], pts[(j+1)%4], colorScalar, 1, 8 );
}
cv::imwrite(QString("test_lines_rects%1.png").arg(tmpVal++).toLocal8Bit().toStdString(), im);
#endif
return out;
}
std::vector<cv::Vec3i> PreloadSource::getThresholds(const cv::Mat& image,
int& bestThresInd)
{
// the last one is a threshold from lower
_histPositions = 3;
_histWindowWidth = 20;
_histMinPercentage = 0.1;
_minColorHeight = 180;
// Calculate histogram
std::vector<cv::MatND> hist(3);
int histSize = 256; // bin size
float range[] = { 0, 255 };
const float *ranges[] = { range };
std::vector<cv::Mat> channels(3);
cv::split(image, channels);
std::vector<cv::Vec3i> out;
std::vector<float> peakheights;
for (int pos = 0; pos <_histPositions; pos++)
{
cv::Vec3i threshs;
peakheights.push_back(0);
for (auto i=0; i<3; i++)
{
cv::calcHist( &channels[i], 1, 0, cv::Mat(), hist[i], 1, &histSize, ranges, true, false );
cv::normalize(hist[i], hist[i], 0, histSize, cv::NORM_MINMAX, -1, cv::Mat() );
//cv::normalize();
threshs[i] = findHistPeak(hist[i],
true,
_histWindowWidth,
_histMinPercentage,
pos);
if (threshs[i] == -1)
{
threshs[i] = 0;
continue;
}
peakheights[pos] += hist[i].at<float>(threshs[i]);
}
out.push_back(threshs);
}
// determine thresh with highest count
float best=0;
bestThresInd = -1;
for (unsigned int i=0; i<peakheights.size(); i++)
{
if (peakheights[i] > best
&& out[i][0] > _minColorHeight
&& out[i][1] > _minColorHeight
&& out[i][2] > _minColorHeight)
{
best = peakheights[i];
bestThresInd = i;
}
}
// find from black threshold
cv::Vec3i thresh;
for (auto i=0; i<3; i++)
{
cv::calcHist( &channels[i], 1, 0, cv::Mat(), hist[i], 1, &histSize, ranges, true, false );
cv::normalize(hist[i], hist[i], 0, histSize, cv::NORM_MINMAX, -1, cv::Mat() );
thresh[i] = findHistPeak(hist[i],
false,
_histWindowWidth,
_histMinPercentage);
if (thresh[i] == -1)
{
thresh[i] = 0;
continue;
}
}
out.push_back(thresh);
return out;
}
std::vector<cv::Vec2f> PreloadSource::extractLines(const std::vector<cv::Vec3i>& thresholds,
const cv::Mat& image)
{
// copy to settings
_imSize = 500;
_histSigma = 15;
_cannyTr1 = 150;
_cannyTr2 = 300;
_medianBlurMask = 15;
_houghRho = 1;
_houghAngle = 180;
_houghFactor = 0.25;
float scale;
cv::Mat im;
if (image.size().width < image.size().height)
{
scale = (float) image.size().height / (float) _imSize;
cv::resize(image, im, cv::Size(int(float(image.size().width)/scale),
_imSize));
}
else
{
scale = (float) image.size().width / (float) _imSize;
cv::resize(image, im, cv::Size(_imSize,
int(float(image.size().height)/scale)));
}
std::vector<cv::Mat> channels(3);
cv::split(im, channels);
// first canny step
cv::Mat canny;
cv::cvtColor(im, canny, cv::COLOR_BGR2GRAY);
cv::Canny(canny, canny, _cannyTr1, _cannyTr2);
/*
* here, different thresholds are tried (meaning that these thresholds are assumed
* as white scanner background) -> more canny steps
*/
cv::Mat mask = cv::Mat::zeros(im.size(), CV_8U);
for (int pos = 0; pos <_histPositions; pos++)
{
// threshold over different histogram peaks (aka pos)
for (auto i=0; i<3; i++)
{
cv::Mat dst(im.size(), CV_8U);
cv::threshold(channels[i], dst, thresholds[pos][i] - _histSigma, 1, cv::THRESH_BINARY_INV);
mask = cv::max(mask, dst);
dst.release();
}
cv::medianBlur(mask, mask, _medianBlurMask);
// apply threshold
std::vector<cv::Mat> newChannels(3);
for (int i=0; i<3; i++)
{
cv::multiply(mask, channels[i], newChannels[i]);
}
cv::Mat tmp;
cv::merge(newChannels, tmp);
cv::cvtColor(tmp, tmp, cv::COLOR_BGR2GRAY);
cv::Mat cannyTmp;
cv::Canny(tmp, cannyTmp, _cannyTr1, _cannyTr2);
canny = cv::max(canny, cannyTmp);
}
/*
* now black scanner background
*/
{
mask = cv::Mat::zeros(im.size(), CV_8U);
for (auto i=0; i<3; i++)
{
cv::Mat dst;
cv::threshold(channels[i], dst, thresholds[_histPositions][i] + _histSigma, 255, cv::THRESH_BINARY_INV);
//
mask = cv::max(mask, dst);
dst.release();
}
cv::medianBlur(mask, mask, _medianBlurMask);
// apply threshold
std::vector<cv::Mat> newChannels(3);
for (int i=0; i<3; i++)
{
newChannels[i] = cv::max(mask, channels[i]);
}
cv::Mat tmp;
cv::merge(newChannels, tmp);
cv::cvtColor(tmp, tmp, cv::COLOR_BGR2GRAY);
cv::Mat cannyTmp;
cv::Canny(tmp, cannyTmp, _cannyTr1, _cannyTr2);
canny = cv::max(canny, cannyTmp );
}
/* continue with hough transform */
//cv::medianBlur(canny, canny, 3); //_medianBlurMaskOverall);
#ifdef DEBUG
cv::imwrite(QString("test_canny%1.png").arg(tmpVal).toLocal8Bit().toStdString(), canny);
#endif
std::vector<cv::Vec2f> lines;
cv::HoughLines(canny, lines, _houghRho, M_PI / float(_houghAngle), int(_houghFactor*_imSize));
for (unsigned int i=0; i<lines.size(); i++)
{
lines[i][0] *= scale;
}
return lines;
}
int PreloadSource::findHistPeak(const cv::Mat& hist,
const bool fromRight,
const int windowWidth,
const float minPercentage,
const int position) const
//hist, fromRight, width = 10, minVal = 0.3, pos = 0):
{
const int n = hist.size[0];
int found = 0;
for (int i = (fromRight ? n-1 : 0); i>=0 && i< n; (fromRight? i-- : i++))
{
bool higher = true;
for (int j=1; j<windowWidth; j++)
{
if (hist.at<float>(i) < hist.at<float>(qMin(i+j, n-1))
|| hist.at<float>(i) < hist.at<float>(qMax(i-j, 0))
|| hist.at<float>(i) < 256.0f * minPercentage)
{
higher = false;
break;
}
}
if (higher)
{
// retrieve pos'th peak
if (found == position)
{
return i;
}
found ++;
}
}
return -1;
}
inline cv::RotatedRect rectFromData(const lbfgsfloatval_t *x)
{
cv::RotatedRect rectangle;
rectangle.angle = x[0];
rectangle.center.x = x[1];
rectangle.center.y = x[2];
rectangle.size.width = x[3];
rectangle.size.height = x[4];
return rectangle;
}
cv::Mat PreloadSource::_currentImageForRectOpt;
lbfgsfloatval_t PreloadSource::evaluate(
void *instance,
const lbfgsfloatval_t *x,
lbfgsfloatval_t *g,
const int n,
const lbfgsfloatval_t stepOrig
)
{
const cv::Mat& image = _currentImageForRectOpt;
cv::RotatedRect rectangle = rectFromData(x);
lbfgsfloatval_t fx = 0.0;
QVector<cv::Point2f> pts(4);
double sumValues;
long int numPixels;
const lbfgsfloatval_t step = stepOrig <= 1e-10 ? 1.0 : stepOrig;
rectangle.points(pts.data());
getSumOfRectangleSamplingF(pts, image, sumValues, numPixels);
fx = -(lbfgsfloatval_t) ( sumValues / (double) numPixels);
// compute gradient
for (int i=0; i<n; i++) {
lbfgsfloatval_t x2[5];
for (int j=0; j < n; j++) {
//x2[j] = j==i ? x[j] + steps[j] : x[j];
x2[j] = j==i ? x[j] + step : x[j];
}
rectangle = rectFromData(x2);
rectangle.points(pts.data());
getSumOfRectangleSamplingF(pts, image, sumValues, numPixels);
lbfgsfloatval_t fx2 = -(lbfgsfloatval_t) ( sumValues / (double) numPixels);
g[i] = (fx2 - fx) / (step );
}
return fx;
}
void PreloadSource::optimizeRectangle(const cv::Mat& edgeMask,
cv::RotatedRect &rectangle)
{
edgeMask.copyTo(_currentImageForRectOpt);
int ret = 0;
const int N = 5;
lbfgsfloatval_t fx;
lbfgsfloatval_t *x = lbfgs_malloc(N);
lbfgs_parameter_t param;
if (x == NULL) {
qCritical() << "ERROR: Failed to allocate a memory block for variables.";
return;
}
/* Initialize the variables. */
x[0] = rectangle.angle;
x[1] = rectangle.center.x;
x[2] = rectangle.center.y;
x[3] = rectangle.size.width;
x[4] = rectangle.size.height;
//qDebug() << "before: " << x[0] << x[1] << x[2] << x[3] << x[4];
/* Initialize the parameters for the L-BFGS optimization. */
lbfgs_parameter_init(&param);
/*param.linesearch = LBFGS_LINESEARCH_BACKTRACKING;*/
param.linesearch = LBFGS_LINESEARCH_BACKTRACKING ;//LBFGS_LINESEARCH_BACKTRACKING_ARMIJO;
/*
Start the L-BFGS optimization; this will invoke the callback functions
evaluate() and progress() when necessary.
*/
ret = lbfgs(N, x, &fx, evaluate, NULL, NULL, &param);
Q_UNUSED(ret);
/* Report the result. */
//qDebug() << "after: " << x[0] << x[1] << x[2] << x[3] << x[4];
//qDebug() << "L-BFGS optimization terminated with status code =" << ret;
rectangle = rectFromData(x);
lbfgs_free(x);
return;
}
Reference in New Issue Copy Permalink