.-- --

スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。
スポンサー広告 comment(-) trackback(-)
.31 2013

ORB特徴点のマッチング

こんにちわ。Uedaです。

OpenCVを使って、ORB特徴点で2枚の画像からどこが同じ特徴のところかを調べてみました。


UIImage *image1 = [UIImage imageNamed:@"IMG_0091.jpg"];
UIImage *image2 = [UIImage imageNamed:@"IMG_0092.jpg"];

Mat imageMat1 = [self cvMatFromUIImage:image1];
Mat imageMat2 = [self cvMatFromUIImage:image2];

Mat grayImage1, grayImage2;
cvtColor(imageMat1, grayImage1, CV_BGR2GRAY);
cvtColor(imageMat2, grayImage2, CV_BGR2GRAY);

OrbFeatureDetector detector(500);
OrbDescriptorExtractor extractor;

std::vector keypoints1, keypoints2;
detector.detect(grayImage1, keypoints1);
detector.detect(grayImage2, keypoints2);

Mat descriptor1, descriptor2;
extractor.compute(grayImage1, keypoints1, descriptor1);
extractor.compute(grayImage2, keypoints2, descriptor2);

vector matches;
BFMatcher matcher(cv::NORM_HAMMING, true);
matcher.match(descriptor1, descriptor2, matches);

for (int i = 0; i < matches.size(); i++)
{
NSLog(@"%f", matches[i].distance);
if (matches[i].distance < 50.0)
{
KeyPoint p1 = keypoints1.at(matches[i].queryIdx);
cv::circle(imageMat1, cv::Point(p1.pt.x, p1.pt.y), 20, Scalar(255, 255, 0), -1);

KeyPoint p2 = keypoints2.at(matches[i].trainIdx);
cv::circle(imageMat2, cv::Point(p2.pt.x, p2.pt.y), 20, Scalar(255, 255, 0), -1);
}
}

[self.imageView1 setImage:[self UIImageFromCVMat:imageMat1]];
[self.imageView2 setImage:[self UIImageFromCVMat:imageMat2]];


結果
orbTest.png

うーん...

以上です。
スポンサーサイト
.30 2013

Eigenでオペレータ演算子を使った行列の計算

こんにちわ。Uedaです。

今回はEigenでオペレータ演算子を使った行列の計算の紹介です。


Matrix3d m1 = Matrix3d::Random(3, 3);
Matrix3d m2 = Matrix3d::Random(3, 3);

Matrix3d m3 = m1 + m2;
Matrix3d m4 = m1 * m2;


直感的に実装できるので便利です。

計算結果
2013-05-30 17:40:33.783 TestEigen[3396:c07]
-0.999984, -0.082700, -0.905911
-0.736924, 0.065534, 0.357729
0.511211, -0.562082, 0.358593
2013-05-30 17:40:33.785 TestEigen[3396:c07]
0.869386, 0.661931, 0.059400
-0.232996, -0.930856, 0.342299
0.038833, -0.893077, -0.984604
2013-05-30 17:40:33.786 TestEigen[3396:c07]
-0.130599, 0.579231, -0.846510
-0.969920, -0.865321, 0.700028
0.550043, -1.455158, -0.626011
2013-05-30 17:40:33.787 TestEigen[3396:c07]
-0.885282, 0.224109, 0.804256
-0.642049, -0.868276, -0.373563
0.589327, 0.541352, -0.515106

以上です。
.29 2013

ORB特徴点を表示してみた

こんにちわ。Uedaです。

前回, iphone側で動画像を取得して, それをOpenCVで扱える形に変換するところまで紹介しました。
今回はとりあえずORBを使って特徴点を表示してみました。


OrbFeatureDetector detector(500);
// OrbDescriptorExtractor *extractor;
std::vector keypoints;
detector.detect(grayImage, keypoints);

for (int i = 0; i < keypoints.size(); i++)
{
KeyPoint p = keypoints.at(i);
cv::circle(imageMat, cv::Point(p.pt.x, p.pt.y), 1, Scalar(0, 255, 0), -1);
}


これでORB特徴点が表示されます。

今回は以上です。
.28 2013

動画像をOpenCVで解析するための準備

こんにちわ。Uedaです。

今回はタイトル通り, ios上で動画像をOpenCVで解析するための準備をしてみました。

http://reiji1020.hatenablog.com/entry/2012/11/18/165804
のサイト曰く, OpenCV2.4.2からFramework形式のものが既に準備されているので, 導入がすごく簡単になっています。

ライブラリにopencv2.frameworkを追加し, 以下のファイルUIViewControllerを実行すると, 動画像をOpenCVで画像処理・解析するためのMatの形式に変換することができます。

ヘッダファイル

@interface TestOpenCVViewController : UIViewController
{
IBOutlet UIImageView *cameraPreviewImageView;
AVCaptureSession* session;
}

@property(strong, nonatomic) AVCaptureSession *session;
@property(strong, nonatomic) UIImageView *cameraPreviewImageView;

@end


ソースファイル

using namespace cv;

@interface TestOpenCVViewController ()

@end

@implementation TestOpenCVViewController

@synthesize session;
@synthesize cameraPreviewImageView;

- (id)initWithNibName:(NSString *)nibNameOrNil bundle:(NSBundle *)nibBundleOrNil
{
self = [super initWithNibName:nibNameOrNil bundle:nibBundleOrNil];
if (self)
{
// Custom initialization
}
return self;
}

- (void)viewDidLoad
{
[super viewDidLoad];
NSLog(@"viewDidLoad");
// Do any additional setup after loading the view from its nib.
// ビデオキャプチャデバイスの取得
AVCaptureDevice* device;
device = [AVCaptureDevice defaultDeviceWithMediaType:AVMediaTypeVideo];
// オートフォーカスモードを解消
if ([device isFocusPointOfInterestSupported] && [device isFocusModeSupported:AVCaptureFocusModeAutoFocus]) {
NSError *error;
if ([device lockForConfiguration:&error]) {
[device setFocusMode:AVCaptureFocusModeAutoFocus];
[device unlockForConfiguration];
}
}
// デバイス入力の取得
AVCaptureDeviceInput* deviceInput;
deviceInput = [AVCaptureDeviceInput deviceInputWithDevice:device error:NULL];
// ビデオデータ出力の作成
NSMutableDictionary* settings;
AVCaptureVideoDataOutput* dataOutput;
settings = [NSMutableDictionary dictionary];
[settings setObject:[NSNumber numberWithInt:kCVPixelFormatType_32BGRA]
forKey:(id)kCVPixelBufferPixelFormatTypeKey];
dataOutput = [[AVCaptureVideoDataOutput alloc] init];
// [dataOutput autorelease];
dataOutput.videoSettings = settings;
[dataOutput setSampleBufferDelegate:self queue:dispatch_get_main_queue()];
// セッションの作成
self.session = [[AVCaptureSession alloc] init];
// 画質の制御
self.session.sessionPreset = AVCaptureSessionPreset640x480;
[self.session addInput:deviceInput];
[self.session addOutput:dataOutput];

// ステータスバーを隠す
[UIApplication sharedApplication].statusBarHidden = YES;
self.view.bounds = [UIScreen mainScreen].bounds;

// カメラの開始
[self.session startRunning];
}

- (void)didReceiveMemoryWarning
{
[super didReceiveMemoryWarning];
// Dispose of any resources that can be recreated.
}

/***** カメラの動作 *****/
-(void)captureOutput:(AVCaptureOutput *)captureOutput didOutputSampleBuffer:(CMSampleBufferRef)sampleBuffer fromConnection:(AVCaptureConnection *)connection
{
// キャプチャしたフレームからCGImageを作成
UIImage *image = [self imageFromSampleBuffer:sampleBuffer];
cv::Mat imageMat = [self cvMatFromUIImage:image];

/* あとはimageMatをOpenCVで煮るなり焼くなり */

image = [self UIImageFromCVMat:imageMat];
cameraPreviewImageView.image = [UIImage imageWithCGImage:image.CGImage scale:image.scale orientation:UIImageOrientationUp];
}

// サンプルバッファのデータからCGImageRefを生成する
- (UIImage *)imageFromSampleBuffer:(CMSampleBufferRef)sampleBuffer
{
CVImageBufferRef imageBuffer = CMSampleBufferGetImageBuffer(sampleBuffer);
// ピクセルバッファのベースアドレスをロックする
CVPixelBufferLockBaseAddress(imageBuffer, 0);
// Get information of the image
uint8_t *baseAddress = (uint8_t *)CVPixelBufferGetBaseAddressOfPlane(imageBuffer, 0);
size_t bytesPerRow = CVPixelBufferGetBytesPerRow(imageBuffer);
size_t width = CVPixelBufferGetWidth(imageBuffer);
size_t height = CVPixelBufferGetHeight(imageBuffer);
// RGBの色空間
CGColorSpaceRef colorSpace = CGColorSpaceCreateDeviceRGB();

CGContextRef newContext = CGBitmapContextCreate(baseAddress,
width,
height,
8,
bytesPerRow,
colorSpace,
kCGBitmapByteOrder32Little | kCGImageAlphaPremultipliedFirst);
CGImageRef cgImage = CGBitmapContextCreateImage(newContext);
CGContextRelease(newContext);
CGColorSpaceRelease(colorSpace);
CVPixelBufferUnlockBaseAddress(imageBuffer, 0);
UIImage *image = [UIImage imageWithCGImage:cgImage scale:1.0 orientation:UIImageOrientationRight];
CGImageRelease(cgImage);
return image;
}

- (cv::Mat)cvMatFromUIImage:(UIImage *)image
{
CGColorSpaceRef colorSpace = CGImageGetColorSpace(image.CGImage);
CGFloat cols = image.size.width;
CGFloat rows = image.size.height;

cv::Mat cvMat(rows, cols, CV_8UC4); // 8 bits per component, 4 channels

CGContextRef contextRef = CGBitmapContextCreate(cvMat.data, // Pointer to data
cols, // Width of bitmap
rows, // Height of bitmap
8, // Bits per component
cvMat.step[0], // Bytes per row
colorSpace, // Colorspace
kCGImageAlphaNoneSkipLast |
kCGBitmapByteOrderDefault); // Bitmap info flags

CGContextDrawImage(contextRef, CGRectMake(0, 0, cols, rows), image.CGImage);
CGContextRelease(contextRef);
// CGColorSpaceRelease(colorSpace);

return cvMat;
}

- (UIImage *)UIImageFromCVMat:(cv::Mat)cvMat
{
NSData *data = [NSData dataWithBytes:cvMat.data length:cvMat.elemSize()*cvMat.total()];
CGColorSpaceRef colorSpace;

if (cvMat.elemSize() == 1) {
colorSpace = CGColorSpaceCreateDeviceGray();
} else {
colorSpace = CGColorSpaceCreateDeviceRGB();
}

CGDataProviderRef provider = CGDataProviderCreateWithCFData((__bridge CFDataRef)data);

// Creating CGImage from cv::Mat
CGImageRef imageRef = CGImageCreate(cvMat.cols, //width
cvMat.rows, //height
8, //bits per component
8 * cvMat.elemSize(), //bits per pixel
cvMat.step[0], //bytesPerRow
colorSpace, //colorspace
kCGImageAlphaNone|kCGBitmapByteOrderDefault,// bitmap info
provider, //CGDataProviderRef
NULL, //decode
false, //should interpolate
kCGRenderingIntentDefault //intent
);


// Getting UIImage from CGImage
UIImage *finalImage = [UIImage imageWithCGImage:imageRef];
CGImageRelease(imageRef);
CGDataProviderRelease(provider);
CGColorSpaceRelease(colorSpace);

return finalImage;
}


あとはMat形式にした画像をOpenCVで用意されている様々な関数を使って煮るなり, 焼くなりすればOpenCVを堪能することができます。

今回は以上です。
.27 2013

Eigenを使ってSVD(特異値分解)

こんにちわ。Uedaです。

Eigenを使ってSVD(特異値分解)をしてみました。

特異値分解 ~wiki
http://ja.wikipedia.org/wiki/特異値分解


Matrix3d P = Matrix3d::Random(3, 3);
JacobiSVD svd(P, Eigen::ComputeFullU | Eigen::ComputeFullV);
MatrixXd S = svd.singularValues();
Matrix3d U = svd.matrixU();
Matrix3d V = svd.matrixV();


コレスキー分解やLU分解も実装されているので, 便利ですね〜。

今回は以上です。
.24 2013

Eigenを使ってQR分解

こんにちわ。Uedaです。

前回, 算出した回転行列の直行行列としての性質が保証されていなかったので, QR分解します。

直行行列 ~wiki
http://ja.wikipedia.org/wiki/直交行列

QR分解は行列を直行行列と上三角行列に分解する行列です。なので, 前回算出した回転行列もどきをQR分解して, 直行行列とそうでない行列に分けて, 直行行列の方を回転行列として抽出しようということです。

回転行列もどきRをQR分解して, 直行行列を抽出します。Eigen使ったコードです。

HouseholderQR QR(R);
Matrix3d Q = QR.householderQ();


出力結果
QR分解前
0.997029, 0.077192, -0.510096
-0.015336, 1.234692, -0.650774
-0.119096, 0.023709, 0.528783

QR分解後
-0.992825, -0.018295, 0.118166
0.015271, -0.999534, -0.026448
0.118594, -0.024454, 0.992642

直行行列にできたはいいですが, 符号がバラバラになってしまってる...
.23 2013

Eigenで最小2乗法を解いてみた 〜再投影誤差編

こんにちわ。Uedaです。

今回は下の資料を参考にカメラ運動パラメータを推定する計算をしてみました。

http://www.ime.info.hiroshima-cu.ac.jp/~hiura/lec/iip/geometry.pdf

最小2乗法の計算はEigenのjacobiSVDを使いました。


const int DATA_COUNT = 100;
const double ANS_P11 = 1.0, ANS_P12 = 0.0, ANS_P13 = 0.0, ANS_P14 = 1.0;
const double ANS_P21 = 0.0, ANS_P22 = 1.0, ANS_P23 = 0.0, ANS_P24 = 1.0;
const double ANS_P31 = 0.0, ANS_P32 = 0.0, ANS_P33 = 1.0, ANS_P34 = 1.0;

const double NOIZE1 = 1.0;

double *Xa = new double[DATA_COUNT];
double *Ya = new double[DATA_COUNT];
double *Za = new double[DATA_COUNT];
double *ua = new double[DATA_COUNT];
double *va = new double[DATA_COUNT];

for (int i = 0; i < DATA_COUNT; i++)
{
double X = 0.5 - [self random];
double Y = 0.5 - [self random];
double Z = 0.5 - [self random];
double u = (ANS_P11*X + ANS_P12*Y + ANS_P13*Z + ANS_P14) / (ANS_P31*X + ANS_P32*Y + ANS_P33*Z + ANS_P34) + (0.5 - [self random])*NOIZE1;
double v = (ANS_P21*X + ANS_P22*Y + ANS_P23*Z + ANS_P24) / (ANS_P31*X + ANS_P32*Y + ANS_P33*Z + ANS_P34) + (0.5 - [self random])*NOIZE1;
Xa[i] = X;
Ya[i] = Y;
Za[i] = Z;
ua[i] = u;
va[i] = v;
}

MatrixXd A(DATA_COUNT*2, 11);
VectorXd b(DATA_COUNT*2);

for (int i = 0, j = 0; j < DATA_COUNT; i+=2, j++)
{
b(i) = ua[j];
b(i+1) = va[j];

A(i, 0) = Xa[j]; A(i, 1) = Ya[j]; A(i, 2) = Za[j]; A(i, 3) = 1.0;
A(i, 4) = 0.0; A(i, 5) = 0.0; A(i, 6) = 0.0; A(i, 7) = 0.0;
A(i, 8) = -Xa[j]*ua[j]; A(i, 9) = -Ya[j]*ua[j]; A(i, 10) = -Za[j]*ua[j];

A(i+1, 0) = 0.0; A(i+1, 1) = 0.0; A(i+1, 2) = 0.0; A(i+1, 3) = 0.0;
A(i+1, 4) = Xa[j]; A(i+1, 5) = Ya[j]; A(i+1, 6) = Za[j]; A(i+1, 7) = 1.0;
A(i+1, 8) = -Xa[j]*va[j]; A(i+1, 9) = -Ya[j]*va[j]; A(i+1, 10) = -Za[j]*va[j];
}

MatrixXd ans = A.jacobiSvd(ComputeThinU | ComputeThinV).solve(b);
// MatrixXd ans = (A.transpose() * A).inverse() * A.transpose() * b;

NSLog(@"\n %f, %f, %f, %f \n %f, %f, %f, %f \n %f, %f, %f"
, ans(0, 0), ans(1, 0), ans(2, 0), ans(3, 0),
ans(4, 0), ans(5, 0), ans(6, 0), ans(7, 0),
ans(8, 0), ans(9, 0), ans(10, 0));


正解データ
1.0 0.0 0.0 1.0
0.0 1.0 0.0 1.0
0.0 0.0 1.0 1.0

出力結果
0.830548 -0.178821 -0.476796 1.108733
-0.103757 0.920693 -0.573100 1.011145
-0.072787 -0.137976 0.506778 1.0

前回計算で失敗した原因は, 関数を係数を定数倍しても意味が不変であることを考慮していなかったためだと思います。今回は参考資料して, そこのところを解消しました。
回転行列の性質を満たしていないのでこのまま使うことはできないですが, Eigenの最小2乗法計算手法jacobiSVDのテストも踏まえて計算してみました。
.22 2013

Eigenを使って非線形最小2乗法 ~Levenberg-Marquardt法 2の修正

こんにちわ。Uedaです。

前回のプログラムに誤りがあったので修正します。

一番はじめの関数定義の部分です。

struct reprojection_functor : Functor
{
reprojection_functor(int inputs, int values, double *X, double *Y, double *Z, double *u, double *v)
: inputs_(inputs), values_(values), X(X), Y(Y), Z(Z), u(u), v(v) {}

double *X;
double *Y;
double *Z;
double *u;
double *v;
int operator()(const VectorXd& b, VectorXd& fvec) const
{
for (int i = 0; i < values_; ++i)
{
fvec[i] = (u[i] - (b[0]*X[i] + b[1]*Y[i] + b[2]*Z[i] + b[3]) / (b[8]*X[i] + b[9]*Y[i] + b[10]*Z[i] + b[11])) * (u[i] - (b[0]*X[i] + b[1]*Y[i] + b[2]*Z[i] + b[3]) / (b[8]*X[i] + b[9]*Y[i] + b[10]*Z[i] + b[11]))
+ (v[i] - (b[4]*X[i] + b[5]*Y[i] + b[6]*Z[i] + b[7]) / (b[8]*X[i] + b[9]*Y[i] + b[10]*Z[i] + b[11])) * (v[i] - (b[4]*X[i] + b[5]*Y[i] + b[6]*Z[i] + b[7]) / (b[8]*X[i] + b[9]*Y[i] + b[10]*Z[i] + b[11]));
}
return 0;
}

const int inputs_;
const int values_;
int inputs() const { return inputs_; }
int values() const { return values_; }
};


なお、前回のテストケースでの実験ではなかなかいい結果を得ることができませんでした。ということで、少しのノイズを軽減させ、正解データも変えました。

const int DATA_COUNT = 100;
const double ANS_P11 = 1.0, ANS_P12 = 0.0, ANS_P13 = 0.0, ANS_P14 = 1.0;
const double ANS_P21 = 0.0, ANS_P22 = 1.0, ANS_P23 = 0.0, ANS_P24 = 1.0;
const double ANS_P31 = 0.0, ANS_P32 = 0.0, ANS_P33 = 1.0, ANS_P34 = 1.0;

const double NOIZE1 = 1.0;
const double NOIZE2 = 1.0;


正しい解に向かっている場合の結果

2013-05-22 17:11:22.906 TestEigen[3914:c07] 正解
1.000000, 0.000000, 0.000000, 1.000000
0.000000, 1.000000, 0.000000, 1.000000
0.000000, 0.000000, 1.000000, 1.000000
2013-05-22 17:11:22.907 TestEigen[3914:c07] 推定前
1.383749, -0.266086, 0.055873, 1.341386
0.277851, 1.302277, 0.388317, 1.318758
-0.484422, -0.412105, 0.862323, 1.235912
2013-05-22 17:11:23.062 TestEigen[3914:c07] 推定後 5
1.048085, -0.024046, 0.006782, 0.932006
0.143338, 0.929105, -0.020849, 0.913913
0.040038, -0.014241, 0.921208, 0.932719


うーん。テストケースをいじっても、3回に1度くらいしか正解らしい結果は得られず、失敗すると大きく誤った結果を返してきます。。

失敗した結果

2013-05-22 17:22:22.372 TestEigen[3956:c07] 正解
1.000000, 0.000000, 0.000000, 1.000000
0.000000, 1.000000, 0.000000, 1.000000
0.000000, 0.000000, 1.000000, 1.000000
2013-05-22 17:22:22.374 TestEigen[3956:c07] 推定前
1.151559, 0.088669, -0.113892, 1.024381
-0.037348, 0.584864, 0.104453, 1.257278
0.182822, 0.103723, 0.785262, 0.762810
2013-05-22 17:22:22.542 TestEigen[3956:c07] 推定後 5
88.457310, -9.044245, -13.995677, 84.093021
1.123463, 78.481974, -10.546093, 81.876606
5.283922, -2.218911, 70.051671, 80.691307
.21 2013

Eigenを使って非線形最小2乗法 ~Levenberg-Marquardt法 2

こんにちわ。Uedaです。

今回は ui = (Xi*p11 + Yi*P12 + Zi*P13 + P14) / (Xi*p31 + Yi*P32 + Zi*P33 + P34)
vi = (Xi*p21 + Yi*P22 + Zi*P23 + P24) / (Xi*p31 + Yi*P32 + Zi*P33 + P34)
の関数についてLevenberg-Marquardt法を使ってみます。

ちなみに、この関数は動画像から再投影誤差を最小にすることで、カメラの運動を推定するときに使ったりします。

前回の関数の定義を

reprojection_functor(int inputs, int values, double *X, double *Y, double *Z, double *u, double *v)
: inputs_(inputs), values_(values), X(X), Y(Y), Z(Z), u(u), v(v) {}

double *X;
double *Y;
double *Z;
double *u;
double *v;
int operator()(const VectorXd& b, VectorXd& fvec) const
{
for (int i = 0; i < values_; ++i)
{
fvec[i] = u[i] - ((b[0]*X[i] + b[1]*Y[i] + b[2]*Z[i] + b[3]) / (b[8]*X[i] + b[9]*Y[i] + b[10]*Z[i] + b[11]))
+ v[i] - ((b[4]*X[i] + b[5]*Y[i] + b[6]*Z[i] + b[7]) / (b[8]*X[i] + b[9]*Y[i] + b[10]*Z[i] + b[11]));
}
return 0;
}

に変更して、テストケースを作れば完成です。

テストケース

const int DATA_COUNT = 100;
const double ANS_P11 = 1.0, ANS_P12 = 4.0, ANS_P13 = 7.0, ANS_P14 = 10.0;
const double ANS_P21 = 2.0, ANS_P22 = 5.0, ANS_P23 = 8.0, ANS_P24 = 11.0;
const double ANS_P31 = 3.0, ANS_P32 = 6.0, ANS_P33 = 9.0, ANS_P34 = 12.0;

const double NOIZE1 = 2.0;
const double NOIZE2 = 2.0;

double *Xa = new double[DATA_COUNT];
double *Ya = new double[DATA_COUNT];
double *Za = new double[DATA_COUNT];
double *ua = new double[DATA_COUNT];
double *va = new double[DATA_COUNT];

double p11, p12, p13, p14, p21, p22, p23, p24, p31, p32, p33, p34;
for (int i = 0; i < DATA_COUNT; i++)
{
double X = 0.5 - [self random];
double Y = 0.5 - [self random];
double Z = 0.5 - [self random];
double u = (ANS_P11*X + ANS_P12*Y + ANS_P13*Z) / (ANS_P31*X + ANS_P32*Y + ANS_P33*Z) + (0.5 - [self random])*NOIZE1;
double v = (ANS_P21*X + ANS_P22*Y + ANS_P23*Z) / (ANS_P31*X + ANS_P32*Y + ANS_P33*Z) + (0.5 - [self random])*NOIZE1;
Xa[i] = X;
Ya[i] = Y;
Za[i] = Z;
ua[i] = u;
va[i] = v;
}

VectorXd p(12);
p11 = ANS_P11 + (0.5 - [self random])*NOIZE2; p12 = ANS_P12 + (0.5 - [self random])*NOIZE2; p13 = ANS_P13 + (0.5 - [self random])*NOIZE2; p14 = ANS_P14 + (0.5 - [self random])*NOIZE2;
p21 = ANS_P21 + (0.5 - [self random])*NOIZE2; p22 = ANS_P22 + (0.5 - [self random])*NOIZE2; p23 = ANS_P23 + (0.5 - [self random])*NOIZE2; p24 = ANS_P24 + (0.5 - [self random])*NOIZE2;
p31 = ANS_P31 + (0.5 - [self random])*NOIZE2; p32 = ANS_P32 + (0.5 - [self random])*NOIZE2; p33 = ANS_P33 + (0.5 - [self random])*NOIZE2; p34 = ANS_P34 + (0.5 - [self random])*NOIZE2;
p(0) = p11; p(1) = p12; p(2) = p13; p(3) = p14;
p(4) = p21; p(5) = p22; p(6) = p23; p(7) = p24;
p(8) = p31; p(9) = p32; p(10) = p33; p(11) = p34;


結果
正解
1.000000, 4.000000, 7.000000, 10.000000
2.000000, 5.000000, 8.000000, 11.000000
3.000000, 6.000000, 9.000000, 12.000000
推定
1.001588, 4.518335, 6.761741, 10.181617
2.237630, 4.923387, 7.199742, 10.156541
3.348539, 6.064186, 9.443100, 12.882333

以上です。
.20 2013

Eigenを使って非線形最小2乗法 ~Levenberg-Marquardt法 補足

こんにちわ。Uedaです。

前回の補足します。

<問題> z = a / x + b / y + c の非線形の式のa, b, cを、ノイズを乗せたM個のデータxi, yi, zi(0 < i < M)から推定する。
・正解a, b, cはそれぞれ
a = 1.0 ,b = 2.0, c = 3.0
・初期データ_a, _b, _cはそれぞれ
_a = a + noise2, _b = b + noise2, _b = b + noise
・-0.5 < xi < 0.5, -0.5 < yi < 0.5, zi = z + noise1
・-2.5 < noise1 < 2.5
・-50.0 < noise2 < 50.0

をEigenで実装されているLM法を使って推定します。

まず、データを作っていきます。これらは直接計算に関係ありません。計算の準備です。

const int DATA_COUNT = 100;
const double ANS_A = 1.0;
const double ANS_B = 2.0;
const double ANS_C = 3.0;
const double NOIZE1 = 5.0;
const double NOIZE2 = 100.0;

double *xa = new double[DATA_COUNT];
double *ya = new double[DATA_COUNT];
double *za = new double[DATA_COUNT];
for (int i = 0; i < DATA_COUNT; i++)
{
double x = 0.5 - [self random];
double y = 0.5 - [self random];
double noize = (0.5 - [self random]) * NOIZE1;
double z = ANS_A/x + ANS_B/y + ANS_C + noize;
xa[i] = x;
ya[i] = y;
za[i] = z;
}


初期値を適当に決定します。今回は正解データにnoise2を足します。(ただし, -50.0 < noise2 < 50.0

VectorXd p(3); // beta1とbeta2の初期値(適当)
double a = ANS_A + (0.5 - [self random]) * NOIZE2;
double b = ANS_B + (0.5 - [self random]) * NOIZE2;
double c = ANS_C + (0.5 - [self random]) * NOIZE2;
p(0) = a;
p(1) = b;
p(2) = c;


関数を定義します。
int operator()(const VectorXd& b, VectorXd& fvec) constの第一引数 bは求めたい係数です。今回ではa, b, cにあたります。初期計算では初期値がこれになります。fvec[i]はあるデータi(xi, yi, zi)に対して, 推定中の関数だとどれだけの誤差があるかを入れるものです。

struct misra1a_functor : Functor
{
・・・
int operator()(const VectorXd& b, VectorXd& fvec) const
{
for (int i = 0; i < values_; ++i) {
fvec[i] = z[i] - b[0]/x[i] - b[1]/y[i] - b[2]; // fvec[i] = zi - a / xi - b / yi - c
}
return 0;
}
・・・
}


実際に計算します。

misra1a_functor functor(3, DATA_COUNT, xa, ya, za);

NumericalDiff<misra1a_functor> numDiff(functor);
LevenbergMarquardt<NumericalDiff<misra1a_functor>> lm(numDiff); // 注意 "<"が全角
int info = lm.minimize(p);


1行目: 関数を生成します。
2行目: 生成した関数の(遍)微分した形の関数を用意させています。手で関数の微分をしてもいいのですが, 数値微分させて代用しています。
3, 4行目: Eigenのレーベンバーグ・マーカート法を使って計算させる。結果を配列pに出力しています。

以上です。
.17 2013

Eigenを使って非線形最小2乗法 ~Levenberg-Marquardt法

こんにちわ。Uedaです。

今回は非線形最小二乗法の問題でよく使われるLevenberg-Marquardt法(以下LM法)をEigenを使って解いてみます。
(LM法の理解は筆者は間違っているかもしれません。以下の語尾には「だと思われます」が付きます

参考サイト
http://d.hatena.ne.jp/wildpie/20120712

非線形問題の解き方についてはニュートン法やガウス・ニュートン法があります。
ニュートン法 ~ wiki
http://ja.wikipedia.org/wiki/ニュートン法

ガウス・ニュートン法参考サイト
http://homepage3.nifty.com/skomo/f6/hp6_9.htm

これらの方法は解のにおいのする方向へ探索していき、一歩ずつ解に近づいていくという方法です。(この解のにおいのする方向を調べるのに微分使ってます。

これらの方法は解への収束速度はとても速いのですが、(初期値の取り方を誤ったりすると)誤った解に収束してしまう恐れの大きい方法です。Levenberg-Marquardt法は少し速度を犠牲にする代わりに、この危険性を下げる方法です。
LM法も解のにおいのする方向へ探索していき、一歩ずつ解に近づいていくという方法ですが、その歩幅を工夫しています。LM法の歩幅は最初は小さく、解の方向でありそうだとわかれば次第に歩幅を大きくしていくという具合です。(ガウス・ニュートン法は最初から大股な感じです。

参考サイトをもとにEigenのLM法を試してみました。
<問題> z = a / x + b / y + c の非線形の式のa, b, cを、ノイズを乗せたM個のデータxi, yi, zi(0 < i < M)から推定する。
・正解a, b, cはそれぞれ
a = 1.0 ,b = 2.0, c = 3.0
・初期データ_a, _b, _cはそれぞれ
_a = a + noise2, _b = b + noise2, _b = b + noise
・-0.5 < xi < 0.5, -0.5 < yi < 0.5, zi = z + noise1
・-2.5 < noise1 < 2.5
・-50.0 < noise2 < 50.0

コード



#import "Eigen/Core"
#import "Eigen/Dense"
#import "unsupported/Eigen/NonLinearOptimization"
#import "unsupported/Eigen/NumericalDiff"

・・・

const int DATA_COUNT = 100;
const double ANS_A = 1.0;
const double ANS_B = 2.0;
const double ANS_C = 3.0;
const double NOIZE1 = 5.0;
const double NOIZE2 = 100.0;

double *xa = new double[DATA_COUNT];
double *ya = new double[DATA_COUNT];
double *za = new double[DATA_COUNT];
for (int i = 0; i < DATA_COUNT; i++)
{
double x = 0.5 - [self random];
double y = 0.5 - [self random];
double noize = (0.5 - [self random]) * NOIZE1;
double z = ANS_A/x + ANS_B/y + ANS_C + noize;
xa[i] = x;
ya[i] = y;
za[i] = z;
}

VectorXd p(3); // beta1とbeta2の初期値(適当)
double a = ANS_A + (0.5 - [self random]) * NOIZE2;
double b = ANS_B + (0.5 - [self random]) * NOIZE2;
double c = ANS_C + (0.5 - [self random]) * NOIZE2;
p(0) = a;
p(1) = b;
p(2) = c;

misra1a_functor functor(3, DATA_COUNT, xa, ya, za);

NumericalDiff<misra1a_functor> numDiff(functor);
LevenbergMarquardt<NumericalDiff<misra1a_functor>> lm(numDiff); // 注意 "<"が全角
int info = lm.minimize(p);

NSLog(@"正解 %f, %f, %f", ANS_A, ANS_B, ANS_C);
NSLog(@"推定前 %f, %f, %f", a, b, c);
NSLog(@"推定後 %f, %f, %f, %d", p[0], p[1], p[2], info);
}

template
struct Functor
{
typedef _Scalar Scalar;
enum {
InputsAtCompileTime = NX,
ValuesAtCompileTime = NY
};
typedef Matrix InputType;
typedef Matrix ValueType;
typedef Matrix JacobianType;
};

struct misra1a_functor : Functor
{
misra1a_functor(int inputs, int values, double *x, double *y, double *z)
: inputs_(inputs), values_(values), x(x), y(y), z(z) {}

double *x;
double *y;
double *z;
int operator()(const VectorXd& b, VectorXd& fvec) const
{
for (int i = 0; i < values_; ++i) {
fvec[i] = z[i] - b[0]/x[i] - b[1]/y[i] - b[2];
}
return 0;
}

const int inputs_;
const int values_;
int inputs() const { return inputs_; }
int values() const { return values_; }
};



結果

2013-05-17 17:28:55.085 TestEigen[2043:c07] 正解 1.000000, 2.000000, 3.000000
2013-05-17 17:28:55.086 TestEigen[2043:c07] 推定前 44.168610, -46.338116, 0.831853
2013-05-17 17:28:55.086 TestEigen[2043:c07] 推定後 0.996553, 1.992320, 2.897285, 1

むちゃくちゃな初期値(結果の推定前)から、見事、正解に近づくことができました!
次回コードレビューします。
.15 2013

Eigenで最小2乗法(3) 〜ロバスト推定(補足)

こんにちわ。Uedaです。

前回は最小2乗法のロバスト性を高める方法、M推定のTukeyのBiweight推定法を紹介しました。
TukeyのBiweight推定法は重みの付け方が重要となりますが、前回はここをはしょってしまいましたので、今回はそこをすこしだけ説明します。
簡単のため、y = ax+bのa, bを推定する問題を考えます。 誤差を含むデータa, bと観測データxi, yi(0 < i < M)からa, bを推定します。
このとき、TukeyのBiweight推定法でも重みはi番目の観測データにたいし、以下で算出されるwとなります。またWは事前に用意しておくもので、これ以上推定している直線から離れている観測データは使いたくありませんっていう値です。
・d = y - ax - b
・w = (1 - (d / W)^2)^2
です。
この重みを採用することで、誤差にひきずられにくい推定ができます。
.14 2013

Eigenで最小2乗法(2) 〜ロバスト推定

こんにちわ。Uedaです。

最小2乗法について、最初の方で触れました。前回、たくさんの観測データ(ノイズを含む)から、その観測データの傾向を抽出するという最適化問題として、利用しました。

今回は観測データに大きめのノイズがのったデータを含ませます。このノイズだらけのデータから、ノイズの影響をあまり受けないように推定する方法を検証します。

ちなみに、ノイズの影響をあまり受けずに推定する、ノイズに強い推定のことを、ロバストな推定、ロバスト推定といいます。

扱うのはTukeyのBiweight推定法です。
TukeyのBiweight推定法について
http://imagingsolution.blog107.fc2.com/blog-entry-32.html

では、検証実験内容です。
w=ax+by+cz+dの式の係数a, b, c, dをノイズの乗ったデータw, x, y, zから推定する。(a, b, c, dは1.0, 2.0, 3.0, 4.0とした) データw, x, y, zはそれぞれ130個あり、
100個は-0.5〜0.5のノイズをのせる。
30個は-50〜50のノイズをのせる。

TukeyのBiweight推定法は重み付け最小2乗法を繰り返し行って、解に近づけていきます。推定した係数をもとに、あるデータがその係数からなる線からどれだけ遠いか(どれだけの誤差があるか)を調べ、近いデータ(誤差が小さい)ほど重みを軽くし、次の最小2乗法での計算で重視します。(筆者の理解が正しければ...
今回繰り返し回数は3回にし、重みを決定するための値を40.0としました。

以下コードと出力です。


const double ANS_A = 1.0;
const double ANS_B = 2.0;
const double ANS_C = 3.0;
const double ANS_D = 4.0;
const int DATA_COUNT_A = 100;
const int DATA_COUNT_B = 30;
const int DATA_COUNT = DATA_COUNT_A + DATA_COUNT_B;
const double NOISE_A = 1.0;
const double NOISE_B = 100.0;
const int CALC_ITERATE_COUNT = 3;
const double W = 40.0;

double* xArr = new double[DATA_COUNT];
double* yArr = new double[DATA_COUNT];
double* zArr = new double[DATA_COUNT];
double* wArr = new double[DATA_COUNT];

double* weightArr = new double[DATA_COUNT];

for (int i = 0; i < DATA_COUNT_A; i++) {
double x = 0.5 - [self random];
double y = 0.5 - [self random];
double z = 0.5 - [self random];
double noize = (0.5 - [self random]) * NOISE_A;
double w = x * ANS_A + y * ANS_B + z * ANS_C + ANS_D + noize;
wArr[i] = w;
xArr[i] = x;
yArr[i] = y;
zArr[i] = z;
weightArr[i] = 1.0;
NSLog(@"%f, %f, %f, %f, %f", x, y, z, w, noize);
}

for (int i = DATA_COUNT_A; i < DATA_COUNT; i++) {
double x = 0.5 - [self random];
double y = 0.5 - [self random];
double z = 0.5 - [self random];
double noize = (0.5 - [self random]) * NOISE_B;
double w = x * ANS_A + y * ANS_B + z * ANS_C + ANS_D + noize;
wArr[i] = w;
xArr[i] = x;
yArr[i] = y;
zArr[i] = z;
weightArr[i] = 1.0;
NSLog(@"%f, %f, %f, %f, %f", x, y, z, w, noize);
}


MatrixXd X(DATA_COUNT, 4);
MatrixXd Y(DATA_COUNT, 1);
MatrixXd ans;

for (int j = 0; j < CALC_ITERATE_COUNT; j++)
{
for (int i = 0; i < DATA_COUNT; i++)
{
X(i, 0) = xArr[i] * weightArr[i];
X(i, 1) = yArr[i] * weightArr[i];
X(i, 2) = zArr[i] * weightArr[i];
X(i, 3) = 1.0 * weightArr[i];
Y(i, 0) = wArr[i] * weightArr[i];
}

MatrixXd Xt = X.transpose();
ans = (Xt * X).inverse() * Xt * Y;

for (int i = 0; i < DATA_COUNT; i++)
{
double ax = xArr[i] * ans(0, 0);
double by = yArr[i] * ans(1, 0);
double cz = zArr[i] * ans(2, 0);
double d = ans(3, 0);
double w = wArr[i];
double e = w - ax - by - cz - d;
// NSLog(@"%f", e);
if (fabs(e) > W) weightArr[i] = 0.0;
else
{
e /= W;
weightArr[i] = (1.0 - e*e) * (1.0 - e*e);
}
}
NSLog(@"推定 %f, %f, %f, %f", ans(0, 0), ans(1, 0), ans(2, 0), ans(3, 0));
}

NSLog(@"正解 %f, %f, %f, %f", ANS_A, ANS_B, ANS_C, ANS_D);


出力
Msuitei2.png


一番目の推定結果はふつうに最小2乗法を解いた場合です。繰り返し計算することにより解に近づいていることがわかります。


.13 2013

画像処理 〜HarrisとStephens、Plesseyのコーナー検出アルゴリズム 2

こんにちわ。Uedaです。

前回特徴点抽出の実験として、HarrisとStephens、Plesseyのコーナー検出アルゴリズムを実装し検証してみました。今回はアルゴリズムの内容をみていこうかと思います。

ただし、筆者の理解は間違ってる可能性があるので、その点ご了承ください。

wiki
http://ja.wikipedia.org/wiki/コーナー検出法

0. 対象とする画像をIとする。Iをグレースケール化する。
1. 大きさ(u, v)のパッチPを用意する(一般的に考えてパッチPは画像Iより(かなり)小さい)。
2. Iの画像上のある位置にパッチPを置く。このパッチの下になっているIの画素をI_P1とする。
3. 次にパッチを(x, y)だけずらした位置に置く(一般的にわずかにずらす)。 このパッチの下になっているIの画素をI_P2とする。
4. I_P1とIP_2の二つを比較して、そこに特徴点が存在するかどうか調べる。
調べ方はI_P1とI_P2の画素の勾配に着目する。(具体的には、Iはグレースケールなので、黒→白となっていれば勾配は大きいとなる )。ここで勾配をもとに作ったものがwikiに記載されている構造テンソルAである(と思われる。
5. Aの固有値を調べ大きかったら特徴点があると判断。(固有値の大きさを使うやり方は主成分分析などでもよく使われています。

2. 〜 5. をIの全画素上で行う。
という手順です。(と思われます。

以上です!!
どうですか?実装できそうではないでしょうか?
.10 2013

画像処理 〜HarrisとStephens、Plesseyのコーナー検出アルゴリズム

こんにちわ。Uedaです。

今日は誰しもOpenCVで一度はお世話になったことのある、コーナー(特徴点)検出についてです。

wiki
http://ja.wikipedia.org/wiki/コーナー検出法

OpenCVでは「goodFeaturesToTrack」関数の前身にあたると思います。筆者は大学時代から今に至るまで、この関数にお世話なりっぱなしです。
ということで、こいつはいったいなにをしているのか、特徴点抽出とはいったい何をやっているのかを簡単に調べてみました。

参考
http://d.hatena.ne.jp/michellon/20120304/1330850172


コード

前回コードで白黒画像を作って...


// 微分画像の生成
std::vector> differential_x;
std::vector> differential_y;
for (y = 1; y < height - 1; y++) {
std::vector dx;
std::vector dy;
for (x = 1; x < width - 1; x++) {

UInt8 *g11, *g12, *g13, *g21, *g22, *g23, *g31, *g32, *g33;
g11 = buffer + (y-1) * bytesPerRow + (x-1)*4; g12 = buffer + (y-1) * bytesPerRow + (x )*4; g13 = buffer + (y-1) * bytesPerRow + (x+1)*4;
g21 = buffer + (y ) * bytesPerRow + (x-1)*4; g22 = buffer + (y ) * bytesPerRow + (x )*4; g23 = buffer + (y ) * bytesPerRow + (x+1)*4;
g31 = buffer + (y+1) * bytesPerRow + (x-1)*4; g32 = buffer + (y+1) * bytesPerRow + (x )*4; g33 = buffer + (y+1) * bytesPerRow + (x+1)*4;

// NSLog(@"%d, %d", *g21 - *g23, *g12 - *g32);
dx.push_back(*g21 - *g23);
dy.push_back(*g12 - *g32);
}
differential_x.push_back(dx);
differential_y.push_back(dy);
}

// 各画素において特徴点かどうか判定
std::vector corners;
for (y = 3; y < height - 3; y++)
{
for (x = 3; x < width - 3; x++)
{
UInt8 A = 0, B = 0, C = 0;
UInt8 x11, x12, x13, x21, x22, x23, x31, x32, x33;
UInt8 y11, y12, y13, y21, y22, y23, y31, y32, y33;

x11 = (differential_x.at(y-1)).at(x-1); x12 = (differential_x.at(y-1)).at( x ); x13 = (differential_x.at(y-1)).at(x+1);
x21 = (differential_x.at( y )).at(x-1); x22 = (differential_x.at( y )).at( x ); x23 = (differential_x.at( y )).at(x+1);
x31 = (differential_x.at(y+1)).at(x-1); x32 = (differential_x.at(y+1)).at( x ); x33 = (differential_x.at(y+1)).at(x+1);

y11 = (differential_y.at(y-1)).at(x-1); y12 = (differential_y.at(y-1)).at( x ); y13 = (differential_y.at(y-1)).at(x+1);
y21 = (differential_y.at( y )).at(x-1); y22 = (differential_y.at( y )).at( x ); y23 = (differential_y.at( y )).at(x+1);
y31 = (differential_y.at(y+1)).at(x-1); y32 = (differential_y.at(y+1)).at( x ); y33 = (differential_y.at(y+1)).at(x+1);

A = x11*x11 + x12*x12 + x13*x13
+ x21*x21 + x22*x22 + x23*x23
+ x31*x31 + x32*x32 + x33*x33;

B = y11*y11 + y12*y12 + y13*y13
+ y21*y21 + y22*y22 + y23*y23
+ y31*y31 + y32*y32 + y33*y33;

C = x11*y11 + x12*y12 + x13*y13
+ x21*y21 + x22*y22 + x23*y23
+ x31*y31 + x32*y32 + x33*y33;

int det = A*B - C*C;
int tr = A + B;

// 判定 (0.1←感度k(sensitivity) 0.04〜0.15
if (35000 < det - 0.1 * tr * tr)
{
corners.push_back(Vector2d(x, y));
}
}
}

// 表示
for (int i = 0; i < corners.size(); i++)
{
Vector2d vec = corners.at(i);
NSLog(@"%f, %f", vec(0), vec(1));
}


検証データ

harris.png

そういえば、感度kってgoodFeaturesToTrackの引数にあったなー!!ってことくらいしか理解できてませんw
.09 2013

画像処理 〜Prewittでエッジ抽出

こんにちわ。Uedaです。

前回Sobelフィルタを使ってエッジ抽出しましたが、今回はPrewittフィルタを使います。

Sobelと基本的に同じだと思います。Sobelでは勾配の検出するためのフィルタに重みを付けたものですが、これをなくしたものだと思ってます。

ということで、Sobelと違うとこだけ、コードあげときます。


// Sobel
// UInt8 gHs = - *g11 + *g13
// - *g21 * 2.0 + *g23 * 2.0
// - *g31 + *g33;
// UInt8 gVs = - *g11 - *g12 * 2.0 - *g13
// + *g31 + *g32 * 2.0 + *g33;

// Prewitt
UInt8 gHs = - *g11 + *g13
- *g21 + *g23
- *g31 + *g33;
UInt8 gVs = - *g11 - *g12 - *g13
+ *g31 + *g32 + *g33;


検証データ
Prewitt.png

前回のSobelフィルタの結果と、あまり違いがありませんねw Sobelの方がくっきりしたエッジになるのではと思っていたのですが...
.08 2013

画像処理 〜Sobelでエッジ抽出

こんにちわ。Uedaです。

前回、前々回で、画像白黒したり、画像をぼかしながら画素を直接いじったりしてみたのは、エッジ抽出や、特徴点抽出の勉強をしたかったためです!

今回はSobelを使ったエッジ抽出をやってみました。

参考サイト
http://www.mis.med.akita-u.ac.jp/~kata/image/sobelprew.html

ある画素Piの周りの勾配を調べ、勾配が大きかったらPiはエッジだねってことで白くしましょってな感じで理解。

コード


std::vector> pixels;
for (y = 1; y < height - 1; y++) {
std::vector rowpixel;
for (x = 1; x < width - 1; x++) {

///////////////////////////////////////////////////////////////////////////////////
// ある画素P(x, y)の周りの画素を取得

UInt8 *g11, *g12, *g13, *g21, *g22, *g23, *g31, *g32, *g33;
g11 = buffer + (y-1) * bytesPerRow + (x-1)*4; g12 = buffer + (y-1) * bytesPerRow + (x )*4; g13 = buffer + (y-1) * bytesPerRow + (x+1)*4;
g21 = buffer + (y ) * bytesPerRow + (x-1)*4; g22 = buffer + (y ) * bytesPerRow + (x )*4; g23 = buffer + (y ) * bytesPerRow + (x+1)*4;
g31 = buffer + (y+1) * bytesPerRow + (x-1)*4; g32 = buffer + (y+1) * bytesPerRow + (x )*4; g33 = buffer + (y+1) * bytesPerRow + (x+1)*4;

///////////////////////////////////////////////////////////////////////////////////
// Sobelのフィルタ処理

UInt8 gHs = - *g11 + *g13
- *g21 * 2.0 + *g23 * 2.0
- *g31 + *g33;
UInt8 gVs = - *g11 - *g12 * 2.0 - *g13
+ *g31 + *g32 * 2.0 + *g33;

UInt8 g = sqrt(gHs*gHs + gVs*gVs);

// UInt8 brightness = (*g11 + *g12 + *g13 + *g21 + *g22 + *g23 + *g31 + *g32 + *g33) / 9;

///////////////////////////////////////////////////////////////////////////////////
// ある画素P(x, y)の値の決定

rowpixel.push_back(g);
}
pixels.push_back(rowpixel);
}

///////////////////////////////////////////////////////////////////////////////////
// 表示
for (y = 1; y < pixels.size() - 1; y++)
{
std::vector rowpixel = pixels.at(y);
for (x = 1; x < rowpixel.size() - 1; x++)
{
UInt8* tmp;
tmp = buffer + y * bytesPerRow + x * 4;

UInt8 p = rowpixel.at(x);

// NSLog(@"%d", brightness);

*(tmp + 0) = p;
*(tmp + 1) = p;
*(tmp + 2) = p;
}
}


検証

sobel.png


前回使ったビリヤードの画像と見比べて下さい!見事、エッジ抽出できましたー!
(ただ、すんごい遅いです... OpenCVってどうやって高速化してるんだろー
.07 2013

画像処理 ~ぼかし

こんにちわ。Uedaです。

前回、画像処理を行うために画像の2値化を行いました。

今回は「ぼかし」を紹介します。

前回行った2値化のあとに、隣り合う画素で平滑化を行えば完成です。


g11 = buffer + (y-1) * bytesPerRow + (x-1)*4; g12 = buffer + (y-1) * bytesPerRow + (x )*4; g13 = buffer + (y-1) * bytesPerRow + (x+1)*4;
g21 = buffer + (y ) * bytesPerRow + (x-1)*4; g22 = buffer + (y ) * bytesPerRow + (x )*4; g23 = buffer + (y ) * bytesPerRow + (x+1)*4;
g31 = buffer + (y+1) * bytesPerRow + (x-1)*4; g32 = buffer + (y+1) * bytesPerRow + (x )*4; g33 = buffer + (y+1) * bytesPerRow + (x+1)*4;

UInt8 brightness = (*g11 + *g12 + *g13 + *g21 + *g22 + *g23 + *g31 + *g32 + *g33) / 9;



処理後
bokashi.png

前回の画像にぼかし処理を加えました。
ちょっとわかりにくいですが、ぼかしてます。

以上です。
.02 2013

画像処理の下準備 ~画像を白黒に

こんにちわ。Uedaです。

画像処理の下準備として、画像を白黒にしてみました。

参考にしたサイトはこちらになります。
http://blog.livedoor.jp/techblog/archives/65454070.html

白黒にしてみました。
shirokuro.png

次回、画像処理でなにかできればなと思っています。
.01 2013

K-means法理解 (2)

こんにちわ。Uedaです。

今回は、前回のブログの内容の説明をします。

まず、10x10の範囲でランダムな点をDATA_COUNTだけ生成します。これが与えられたデータになります。


for (int i = 0 ; i < DATA_COUNT; i++)
{
double *data = new double[3];
data[0] = [self random] * 10.0;
data[1] = [self random] * 10.0;
double classVal = [self random];
for (int j = 0; j < CLASS_COUNT; j++)
{
if ( classVal < ((j+1) / (double)CLASS_COUNT) )
{
data[2] = j;
break;
}
}
dataArray.push_back(data);
NSLog(@"%f,%f,%d, %f", data[0], data[1], (int)data[2], classVal);
}


このデータをまとまりに分割しよう(今回は5つ)。分割の方法にk-means法を使おうといった具合です。

前回紹介したサイトに、k-means法の手順がわかりやすく掲載されています。これが、前回実装したコードのどこの部分にあたるか説明します。

① 各点にランダムにクラスタを割り当てる (上のコード)
② クラスタの重心を計算する。
③ 点のクラスタを、一番近い重心のクラスタに変更する
④ 変化がなければ終了。変化がある限りは 2. に戻る。


while (!bFinish)
{

for (int i = 0; i < CLASS_COUNT; i++)
{
double *center = centerArray.at(i);
double x = 0.0, y = 0.0;
int count = 0;
for (int j = 0; j < DATA_COUNT; j++)
{
double *data = dataArray.at(j);
if ((int)data[2] == i) // あるクラスタに割り振られたデータかどうか
{
x += data[0]; // 重心計算中
y += data[1]; //
count++;
}
}
center[0] = x / (double)count; // 重心計算
center[1] = y / (double)count; // 重心計算
}

bFinish = true;

for (int i = 0; i < DATA_COUNT; i++)
{
double *data = dataArray.at(i);
double minDist = 0.0;
int min = -1;
int prev = (int) data[2];
for (int j = 0; j < CLASS_COUNT; j++)
{
double* center = centerArray.at(j);
double dx = data[0] - center[0];
double dy = data[1] - center[1];
double dist = sqrt(dx*dx + dy*dy);
if (min < 0 || dist < minDist) // あるデータがどのクラスタの重心に近いか
{
min = j;
minDist = dist;
data[2] = j;
}
}
int cur = (int) data[2];
if (prev != cur) bFinish = false; // ← ④ (ひとつでもクラスタが変われば、もう一度)
}
calcCount++;
}
NSLog(@"計算回数 %d", calcCount);


以上です。

ようは、DATA_COUNT種類の料理を縦軸に甘い順に、横軸に辛い順に並べて、味の近いかどうかで分割しよう。といったことをやってみました。
 HOME 

ブログ内検索

関連リンク

製品情報

最新記事

カテゴリ

プロフィール

neoxneo



NEXT-SYSTEM iOS Developers Blog


  • UTO:
    カナダ版iPhone4Sは、マナーモードでシャッター音がならない…


  • Ehara:
    ...


  • Hayate:
    ...


  • Tasaki:
    Developer登録完了...したのはいいけど


  • Ueda:
    ...



リンク

上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。