Caffe测试：输入任意一张图片抽取任意一层的feature map，借助MemoryDataLayer

基于extract_feature.cpp修改

extract_feature.cpp是根据prototxt进行测试，因此需要将数据转成特定的格式放入datalayer，一般都是lmdb。

但是在实际应用中，将数据先转成lmdb很麻烦，本文意在实现直接读入3D数据，输入到网络中进行测试，输出一个3D特征图，不需格式转换。

本文的输入数据是三维图像，截成二维的slice进行训练，但是在测试的时候期望能实现输入一张三维图像，直接得到测试后的输出三维feature。因此在以下博文的基础上做了修改

https://blog.csdn.net/u011559236/article/details/78516725

关键由三个部分：

一。数据输入

现将数据读入内存，再使用memorydatalayer将数据放入数据层的Blob中。

memorydatalayer是使用cv::mat读取二维图像，再通过以下代码放入blob

	caffe::MemoryDataLayer<NetF> *layer = (caffe::MemoryDataLayer<NetF> *)feature_extraction_net->layers()[0].get();
	layer->AddMatVector(cv_mat, label);
	

在3D图像中，需要改一下AddMatVector，不要改原函数，以防以后还会用这个函数，可以在memory_data_layer.cpp和memory_data_layer.hpp中加上一个新的函数的定义和声明。

template <typename Dtype>
void MemoryDataLayer<Dtype>::AddRawData(float* mat_vector,
	    const vector<int>& labels, int num) {
	//size_t num = mat_vector.size();
	//size_t num = batch_size_;
	std::cout << "nSlice：" << num << endl;
	std::cout << "nchannels：" << channels_ << endl;
	std::cout << "nheight：" << height_ << endl;
	std::cout << "nwidth：" << width_ << endl;

	CHECK(!has_new_data_) <<
		"Can't add mat until current data has been consumed.";
	CHECK_GT(num, 0) << "There is no mat to add";
	CHECK_EQ(num % batch_size_, 0) <<
		"The added data must be a multiple of the batch size.";
	added_data_.Reshape(num, channels_, height_, width_);
	added_label_.Reshape(num, 1, 1, 1);
	// Apply data transformations (mirror, scale, crop...)
	//this->data_transformer_->Transform(mat_vector, &added_data_);// this layer doesn't need transform //useless in this function, turn it off
	//Dtype* transformed_data = transformed_blob->mutable_cpu_data(); 


	Dtype* feature_blob_data = added_data_.mutable_cpu_data();

	for (int i = 0; i < num*channels_*height_*width_; i++){
	 feature_blob_data[i] = mat_vector[i];
	}// put mat_vector into Blob

这里的mat_vector是一个存储图像数据的数组，不是cv::mat矩阵。

我没有继续用cv::mat，原因是三维数据的格式是rawdata，不能使用opencv读取，若是在AddRawData里沿用cv::mat，还需要将数组转成cv::mat，没有必要。

同时，还要在memory_data_layer.hpp中加上一个声明，一句话

  virtual void AddRawData(float* mat_vector,
  	const vector<int>& labels, int num);

这里仅实现把数据放入Data层的Blob中，没有做data_transform的操作，data transform是在数据层外面做的。

 

二。extract_feature.cpp

读取数据》normalization》interpolation》cropping or padding》extract features》save features

#include <string>
#include <vector>

#include "boost/algorithm/string.hpp"
#include "google/protobuf/text_format.h"

#include "caffe/blob.hpp"
#include "caffe/common.hpp"
#include "caffe/net.hpp"
#include "caffe/layer.hpp"
#include "caffe/proto/caffe.pb.h"
#include "caffe/util/db.hpp"
#include "caffe/util/format.hpp"
#include "caffe/util/io.hpp"
#include "caffe/layers/memory_data_layer.hpp"
//#include  "D:/iyuzu/interpolation/Interpolator.h"
#include  "D:/iyuzu/interpolation/TriLinearInterpolator.h"
//#include "caffe/layers/my_image_data_layer.hpp"

#include <direct.h>
#include <iostream>
#include  <io.h>  
#include <cmath>
#include  <stdio.h>  
#include  <stdlib.h>



#define NetF float 
using caffe::Blob;
using caffe::Caffe;
using caffe::Datum;
using caffe::Net;
using caffe::Layer;
using std::string;
using caffe::MemoryDataLayer;
//using caffe::MyImageDataLayer;
namespace db = caffe::db;

template<typename Dtype>
int feature_extraction_pipeline(int argc, char** argv);
static void CheckFile(const std::string& filename1) {
	std::ifstream f(filename1.c_str());
	if (!f.good()) {
		f.close();
		throw std::runtime_error("Could not open file " + filename1);
	}
	f.close();
}


int main(int argc, char** argv) {
	return feature_extraction_pipeline<float>(argc, argv);
	//  return feature_extraction_pipeline<double>(argc, argv);
}

template<typename Dtype>
int feature_extraction_pipeline(int argc, char** argv) {
	clock_t start = clock();
	::google::InitGoogleLogging(argv[0]);
	const int num_required_args = 12;
	if (argc < num_required_args) {
		LOG(ERROR) <<
			"This program takes in a trained network and an input data layer, and then"
			" extract features of the input data produced by the net.\n"
			"Usage: extract_features  input_image_names pretrained_net_param"
			"  feature_extraction_proto_file  extract_feature_blob_name1[,name2,...]"
			"  save_feature_dataset_name1[,name2,...]  input_size(3 int numbers height width slice)"
			"resolution_of_input(3 int numbers height width slice) input_image_path"
			"  [CPU/GPU] [DEVICE_ID=0]\n"
			"Note: you can extract multiple features in one pass by specifying"
			" multiple feature blob names and dataset names separated by ','."
			" The names cannot contain white space characters and the number of blobs"
			" and datasets must be equal.";
		return 1;
	}
	int arg_pos = num_required_args;

	arg_pos = num_required_args;
	if (argc > arg_pos && strcmp(argv[arg_pos], "GPU") == 0) {
		LOG(ERROR) << "Using GPU";
		int device_id = 0;
		if (argc > arg_pos + 1) {
			device_id = atoi(argv[arg_pos + 1]);
			CHECK_GE(device_id, 0);
		}
		LOG(ERROR) << "Using Device_id=" << device_id;
		Caffe::SetDevice(device_id);
		Caffe::set_mode(Caffe::GPU);
	}
	else {
		LOG(ERROR) << "Using CPU";
		Caffe::set_mode(Caffe::CPU);
	}

	arg_pos = 0;  // the name of the executable
	std::string pretrained_binary_proto(argv[++arg_pos]);//1 caffemodel

	std::string feature_extraction_proto(argv[++arg_pos]);//2 prototxt file
	boost::shared_ptr<Net<Dtype> > feature_extraction_net(
		new Net<Dtype>(feature_extraction_proto, caffe::TEST));
	feature_extraction_net->CopyTrainedLayersFrom(pretrained_binary_proto);


	std::string extract_feature_blob_names(argv[++arg_pos]);//3 blob names in prototxt
	std::vector<std::string> blob_names;
	boost::split(blob_names, extract_feature_blob_names, boost::is_any_of(","));

	std::string save_feature_dataset_names(argv[++arg_pos]);//4 save filename
	std::vector<std::string> dataset_names;
	boost::split(dataset_names, save_feature_dataset_names,
		boost::is_any_of(","));
	CHECK_EQ(blob_names.size(), dataset_names.size()) <<
		" the number of blob names and dataset names must be equal";

	size_t num_features = blob_names.size();

	for (size_t i = 0; i < num_features; i++) { //num_features is refer to the number of features to be extracted
		CHECK(feature_extraction_net->has_blob(blob_names[i]))
			<< "Unknown feature blob name " << blob_names[i]
			<< " in the network " << feature_extraction_proto;
	}

	//int num_mini_batches = atoi(argv[++arg_pos]);//5 num_mini_batches
	int batchsize = 1;// atoi(argv[++arg_pos]);


					  //input size and resolution
	int input_wid = atoi(argv[++arg_pos]);//6
	int input_hei = atoi(argv[++arg_pos]);//7
	int input_slice = atoi(argv[++arg_pos]); //8
	int nVolSize = input_wid*input_hei*input_slice;

	int nWid = 128;// atoi(argv[++arg_pos]);  //defalut width
	int nHei = 128;// atoi(argv[++arg_pos]); 
	int nSlice; //


	double resoX = atof(argv[++arg_pos]); //9 resolution of Width
	double resoY = atof(argv[++arg_pos]); //10 resolution of height
	double resoZ = atof(argv[++arg_pos]); //11 resolution of slice
	double output_VolSize = nHei*nWid*resoZ*input_slice;
	int new_slice = resoZ*input_slice;
	double num_mini_batches = resoZ*input_slice;

	//read input image
	FILE * fpWMov = NULL;
	const char* filename1 = argv[++arg_pos]; //12 input filename 
	CheckFile(filename1);

	fopen_s(&fpWMov, filename1, "rb");
	unsigned short* psData = new unsigned short[input_wid*input_hei*input_slice];
	if (fpWMov)
	{
		fread(psData, sizeof(unsigned short), nVolSize, fpWMov);
		fclose(fpWMov);
	}


	// normalization
	LOG(ERROR) << "Data Processing...";
	LOG(ERROR) << "****Input Data Need Normalization...";
	unsigned short dminvalue = psData[0];
	for (int m = 0; m < input_wid*input_hei*input_slice; m++)
	{
		if (dminvalue > psData[m])
		{
			dminvalue = psData[m];
		}
	}

	unsigned short dmaxvalue = psData[0];
	for (int m = 0; m < input_wid*input_hei*input_slice; m++)
	{
		if (dmaxvalue < psData[m])
		{
			dmaxvalue = psData[m];
		}
	}
	float* pData = new float[input_wid*input_hei*input_slice];
	for (int t = 0; t < nVolSize; t++) {
		pData[t] = (float)(psData[t] - dminvalue) / (dmaxvalue - dminvalue);
	}// normalization
	
	
	// interpolation
	CTriLinearInterpolator<float, float>* interpolator = new CTriLinearInterpolator<float, float>;
	interpolator->SetImgData(pData, input_wid, input_hei, input_slice);

	float *pWriteData = NULL;
	
	pWriteData = new float[new_slice*nWid*nHei];
	if (resoX != 1 || resoY != 1 || resoZ != 1) {
		LOG(ERROR) << "****Input Data Need Interpolation...";
		
		int new_height = resoY*input_hei;
		int new_width = resoX*input_wid;

		double scl_slice = (double)input_slice / new_slice;//dCenterx;
		double scl_width = (double)input_wid / new_width;//dCentery;
		double scl_hei = (double)input_hei / new_height;//dCenterz;
		double dx, dy, dz;

		float * pInterpData = NULL;
		//pInterpData = new unsigned short[new_slice*new_hei*new_width];
		pInterpData = new float[new_slice*new_height*new_width];

		if (!pInterpData)
		{
			return -1;
		}
		for (int n = 0; n < new_slice*new_height*new_width; n++)
		{
			pInterpData[n] = 0;
		}
		float value = 0;
		for (int z = 0; z < new_slice; z++)
		{
			dz = scl_slice*z;
			for (int y = 0; y < new_height; y++)
			{
				dy = scl_hei*y;
				for (int x = 0; x < new_width; x++)
				{
					value = 0;
					dx = scl_width*x;
					interpolator->Interpolate(value, dx, dy, dz);//the matrix is pData
					pInterpData[z*new_height*new_width + y*new_width + x] = value;

				}
			}
		}//interpolation


		 //cropping //z axis needn't crop
		

		if (new_height != nHei || new_width != nWid) {
			float * pCropData = NULL;
			pCropData = new float[new_slice*nWid*nHei];

			for (int n = 0; n < new_slice*nWid*nHei; n++)
			{
				pCropData[n] = 0;
			}
			if (new_height > nHei && new_width>nWid) {
				LOG(ERROR) << "****Input Data Need Cropping...";
				
				double h_off = abs(new_height - nHei) / 2;
				double w_off = abs(new_width - nWid) / 2;

				int top_index_n = 0;
				int data_index_n = 0;
				for (int c = 0; c < new_slice; ++c) {
					int top_index_c = (top_index_n + c) * nHei;
					int data_index_c = (data_index_n + c) * new_height + h_off;
					for (int h = 0; h < nHei; ++h) {
						int top_index_h = (top_index_c + h) * nWid;
						int data_index_h = (data_index_c + h) * new_width + w_off;
						for (int w = 0; w < nWid; ++w) {
							pCropData[top_index_h + w] = pInterpData[data_index_h + w];
						}
					}
				}
			}
			else if (new_height < nHei && new_width < nWid) {
				LOG(ERROR) << "****Input Data Need Padding...";
				double h_off = abs(new_height - nHei) / 2;
				double w_off = abs(new_width - nWid) / 2;
				int top_index_n = 0;
				int data_index_n = 0;
				for (int c = 0; c < new_slice; ++c) {
					int top_index_c = (top_index_n + c) * new_height;
					int data_index_c = (data_index_n + c) * nHei + h_off;
					for (int h = 0; h < new_height; ++h) {
						int top_index_h = (top_index_c + h) * new_width;
						int data_index_h = (data_index_c + h) * nWid + w_off;
						for (int w = 0; w < new_width; ++w) {
							pCropData[data_index_h + w] = pInterpData[top_index_h + w];
						}
					}
				}
			}


			memcpy(pWriteData, pCropData, sizeof(float)*new_slice*nHei*nWid);
			delete[] pCropData;
		}

		else {
			memcpy(pWriteData, pInterpData, sizeof(float)*new_slice*nHei*nWid);
			delete[] pInterpData;
		}
	}
	else {
		memcpy(pWriteData, pData, sizeof(float)*resoZ*input_slice*nHei*nWid);
		delete[] pData;
	}


	std::vector<int> label = { 0 };

	caffe::MemoryDataLayer<NetF> *layer = (caffe::MemoryDataLayer<NetF> *)feature_extraction_net->layers()[0].get();
	layer->AddRawData(pWriteData, label, new_slice);


	// extrat my features

	unsigned char*  Vi = new unsigned char[output_VolSize];//save mask
	unsigned short*  seg = new unsigned short[output_VolSize];//save mask
	int idx = 0;
	LOG(ERROR) << "Extracting Features...";
	
	for (int batch_index = 0; batch_index < num_mini_batches; ++batch_index) {//num_mini_batches
		std::vector<caffe::Blob<NetF>*> input_vec;
		feature_extraction_net->Forward(input_vec);

		for (int i = 0; i < num_features; ++i) {

			const boost::shared_ptr<Blob<Dtype> > feature_blob =
				feature_extraction_net->blob_by_name(blob_names[i]);
			int batch_size = feature_blob->num();
			int dim_features = feature_blob->count() / batch_size;
			const Dtype* feature_blob_data;


			//float*  Vi = new float[nVolSize];
			for (int n = 0; n < batch_size; ++n) {

				feature_blob_data = feature_blob->cpu_data() +
					feature_blob->offset(n);


				for (int d = 0; d < dim_features; ++d) {

					if (feature_blob_data[d]>0.5) {
						Vi[idx] = 255;
						seg[idx] = psData[idx];
					}
					else {
						Vi[idx] = 0;
						seg[idx] = 0;
					}

					++idx;
				}
			}
		}
	}
	//make dir if the save path does not exist;
	string path = dataset_names[0].c_str();
	if (0 != access(path.c_str(), 0))
	{
		mkdir(path.c_str());   
	}

	string seg_file = "\\seg_image.raw";
	string save_path;
	save_path = path + seg_file;
	FILE * fpFeatures = NULL;
	fopen_s(&fpFeatures, save_path.c_str(), "wb+");//93-A2-PDW-mag2_new	
	if (fpFeatures)
	{

		fwrite(seg, sizeof(unsigned short), output_VolSize, fpFeatures);
		fclose(fpFeatures);
	}

	string mask_file = "\\mask_image.raw"; 
	string re_path;
	re_path = path + mask_file;

	FILE * fpmask = NULL;
	fopen_s(&fpmask, re_path.c_str(), "wb+");//93-A2-PDW-mag2_new	
	if (fpmask)
	{

		fwrite(Vi, sizeof(unsigned char), output_VolSize, fpmask);
		fclose(fpmask);
	}

	clock_t end = clock();
	double totaltime;
	totaltime = (double)(end - start) / CLOCKS_PER_SEC;
	std::cout << "Total Time:" << totaltime << "s!" << std::endl;
	LOG(ERROR) << "Successfully extracted the features!";
	return 0;
}

写的比较乱，三维》二维slice》输出二维slice的feature》重组成三维图，这个过程比较复杂，好好看循环应该可以理解的。。。

后半部分就是在把输出的feature 按照一定的顺序存在Vi数组中，存成raw格式，就可以在ImageJ中打开了

输出如图，是三维的，现在是第21个slice。

三。调用格式。

因为我不需要将输出存成lmdb，因此删除了相关代码，因此argv中的倒数第二个参数lmdb也不需要了，刚好我可以把它换成输入图像的路径，因此argc的个数是不变的。

下面是我的cmd命令，改一下路径就好了。

D:\caffe-master\Build\x64\Release\extract_features.exe .\unet_shuf_iter_5600.caffemodel .\train_val_upsamp_pred.prototxt score_out .\features.raw 1 60 .\sub-70074_dwi.raw GPU

最后总结一下吧，看起来很简单的一个功能，但是搞了两个多星期，本来一点都不会写c++，这是第一次接触，刚开始的时候看不懂源码，后来看懂了但是写不出来，因为不知道c++的语法。后来在网上找了很多博主po的代码，一点点对照着写。

先是把存成lmdb的部分改成了存成三维rawdata。输入一开始用的还是lmdb格式的，但是在实际应用中，每次都把三维数据转成lmdb存起来太不实用，下定决心改成直接读取rawdata。

好在发现了memorydatalayer，省了不少麻烦，继续一点一点对照着改。刚开始想的办法是将图像的数组转成三维的cv::mat，也编译成功了，但是结果有问题，因为cv::mat的索引问题，不细说了。后来想直接去掉转cv::mat这一步吧，直接将图像的数组放进datalayer的Blob不就好了，所以加上了AddRawData一层。

曾经觉得自己肯定学不会c++，现在也可以磕磕绊绊写一点了，所以慢慢来吧