/*
# mrcImageMeanFreePathCalc : $Revision$  
# $Date$ 
# Created by $Author$
# Usage : mrcImageMeanFreePathCalc
# Attention
#   $Loccker$
#  	$State$ 
#
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>                  
#define GLOBAL_DECLARATION
#include "../inc/config.h"

#undef DEBUG
#include "genUtil.h"
#include "Memory.h"
#include "mrcImage.h"
#include "lmrcImageMultiplying.h"
#include "lmrcImageTranspose.h"

typedef enum lmrcImageMeanFreePathCalcLineMode {
	lmrcImageMeanFreePathCalcLineModeAnd = 0,
	lmrcImageMeanFreePathCalcLineModeOr  = 1, 
	lmrcImageMeanFreePathCalcLineModeGE  = 2  
} lmrcImageMeanFreePathCalcLineMode; 

typedef enum lmrcImageMeanFreePathCalcLineShape {
	lmrcImageMeanFreePathCalcLineShapeRectangle=0,
	lmrcImageMeanFreePathCalcLineShapeDisk=1
} lmrcImageMeanFreePathCalcLineShape; 

typedef struct lmrcImageMeanFreePathCalcInfo {
	int  N;
	int* countWhite;
	int* countBlack;
	int* countWhiteEdge;
	int* countBlackEdge;
	double averageWhite;
	double averageBlack;
	double averageWhiteEdge;
	double averageBlackEdge;
	int sumWhite;
	int sumBlack;
	int sumWhiteEdge;
	int sumBlackEdge;

	//
	mrcImage  white;
	mrcImage  black;
	mrcImage  whiteEdge;
	mrcImage  blackEdge;

	// Control: Input
	mrcImageParaTypeRealCoord Length;
	mrcImageParaTypeRealCoord LineThickness;   
	lmrcImageMeanFreePathCalcLineMode  LineMode;
	lmrcImageMeanFreePathCalcLineShape LineShape;

	// TemporalVariables
	mrcImageParaTypeRealCoord min; 
	mrcImageParaTypeRealCoord max; 
	mrcImageParaTypeRealCoord cur; 
} lmrcImageMeanFreePathCalcInfo;

typedef enum lmrcImageMeanFreePathCalcMode {
	MeanFreePathCalcModeAlongZaxis=0,
	MeanFreePathCalcModeAlongXaxis=1,
	MeanFreePathCalcModeAlongYaxis=2,
	MeanFreePathCalcModeCalcMode=0x0f,
	MeanFreePathCalcModeRealLength=0x10
} lmrcImageMeanFreePathCalcMode;

extern void lmrcImageMeanFreePathCalcPrint(FILE* fpt, lmrcImageMeanFreePathCalcInfo* linfo, lmrcImageMeanFreePathCalcMode mode);
extern void lmrcImageMeanFreePathCalc(mrcImage* in, lmrcImageMeanFreePathCalcInfo* linfo, lmrcImageMeanFreePathCalcMode mode);
extern void lmrcImageMeanFreePathCalcAlongZaxis(mrcImage* in, lmrcImageMeanFreePathCalcInfo* linfo, lmrcImageMeanFreePathCalcMode  mode);
extern double lmrcImageMeanFreePathCalcBWEvaluation(mrcImage* in, lmrcImageMeanFreePathCalcInfo* linfo, lmrcImageMeanFreePathCalcMode  mode);

extern void lmrcImageMeanFreePathCalcModePrint(FILE* fpt);

int
main(int argc, char* argv[]) 
{
	mrcImageMeanFreePathCalcInfo info;
	lmrcImageMeanFreePathCalcInfo linfo;
	mrcImage in;

	init0(&info);
    argCheck(&info, argc, argv);
    init1(&info);

	DEBUGPRINT("Program Start\n");

	mrcFileRead(&in, info.In, "in main", 0);

	linfo.Length = in.HeaderLength;
	linfo.LineThickness.x = info.LineThicknessX;
	linfo.LineThickness.y = info.LineThicknessY;
	linfo.LineThickness.z = 1;

	linfo.LineMode  = info.LineMode;
	linfo.LineShape = info.LineShape;

	DEBUGPRINT1("mode: %ld\n", info.mode);
	lmrcImageMeanFreePathCalc(&in, &linfo, info.mode);

	lmrcImageMeanFreePathCalcPrint(info.fptOut, &linfo, info.mode);

	if(info.flagOutWhite) {
		mrcFileWrite(&linfo.white, info.OutWhite, "in main", 0);
	}
	if(info.flagOutBlack) {
		mrcFileWrite(&linfo.black, info.OutBlack, "in main", 0);
	}
	if(info.flagOutWhiteEdge) {
		mrcFileWrite(&linfo.whiteEdge, info.OutWhiteEdge, "in main", 0);
	}
	if(info.flagOutBlackEdge) {
		mrcFileWrite(&linfo.blackEdge, info.OutBlackEdge, "in main", 0);
	}
	return EXIT_SUCCESS;
}

void
additionalUsage()
{
	fprintf(stderr, "----- Additional Usage -----\n");
	lmrcImageMeanFreePathCalcModePrint(stderr);
}


void
lmrcImageMeanFreePathCalcModePrint(FILE* fpt)
{
	fprintf(fpt, ">>>> mode \n");
	fprintf(fpt, "%d: AlongZ-axis\n", MeanFreePathCalcModeAlongZaxis);
	fprintf(fpt, "%d: AlongX-axis\n", MeanFreePathCalcModeAlongXaxis);
	fprintf(fpt, "%d: AlongY-axis\n", MeanFreePathCalcModeAlongYaxis);
	fprintf(fpt, "+%d: Length is shown in real using HeaderLength \n", MeanFreePathCalcModeRealLength);

	fprintf(fpt, ">>>> line shape\n");
	fprintf(fpt, "%d: LineShape-Rectangle\n", lmrcImageMeanFreePathCalcLineShapeRectangle);
	fprintf(fpt, "%d: LineShape-Disk\n", lmrcImageMeanFreePathCalcLineShapeDisk);
	fprintf(fpt, ">>>> line mode\n");
	fprintf(fpt, "%d: LineMode-And \n", lmrcImageMeanFreePathCalcLineModeAnd);
	fprintf(fpt, "%d: LineMode-Or  \n", lmrcImageMeanFreePathCalcLineModeOr);
	fprintf(fpt, "%d: LineMode-GE(Greater than or equal to) \n", lmrcImageMeanFreePathCalcLineModeGE);
}

void
lmrcImageMeanFreePathCalc(mrcImage* in, lmrcImageMeanFreePathCalcInfo* linfo, lmrcImageMeanFreePathCalcMode mode)
{
	mrcImage tmp;
	DEBUGPRINT1("mode: %d in Calc 1st \n", mode);
	switch(MeanFreePathCalcModeCalcMode&mode) {
		case MeanFreePathCalcModeAlongZaxis: {
			DEBUGPRINT1("mode: %d in AlongZaxis\n", mode);
			lmrcImageMeanFreePathCalcAlongZaxis(in, linfo, mode);
			break;
		}
		case MeanFreePathCalcModeAlongXaxis: {
			DEBUGPRINT1("mode: %d in AlongXaxis\n", mode);
			lmrcImageTranspose(&tmp, in, 4, 0); 
			lmrcImageMeanFreePathCalcAlongZaxis(&tmp, linfo, mode);
			break;
		}
		case MeanFreePathCalcModeAlongYaxis: {
			DEBUGPRINT1("mode: %d in AlongYaxis\n", mode);
			lmrcImageTranspose(&tmp, in, 2, 0); 
			lmrcImageMeanFreePathCalcAlongZaxis(&tmp, linfo, mode);
			break;
		}
		default: {
			DEBUGPRINT1("mode: %d in default\n", mode);
			fprintf(stderr, "Not supported mode: %d\n", mode);
			fprintf(stderr, "Thus, process as -m %d \n", MeanFreePathCalcModeAlongZaxis);
			lmrcImageMeanFreePathCalcAlongZaxis(in, linfo, mode);
			break;
		}
	}

}

void
lmrcImageMeanFreePathCalcPrint(FILE* fpt, lmrcImageMeanFreePathCalcInfo* linfo, lmrcImageMeanFreePathCalcMode mode)
{
	int i;

	fprintf(fpt, "ANALYSIS Average HowManyLines HowManyVoxels\n");
	fprintf(fpt, "White: %15.6f %d %d\n", linfo->averageWhite, linfo->sumWhite, linfo->countWhite[0]);
	fprintf(fpt, "Black: %15.6f %d %d\n", linfo->averageBlack, linfo->sumBlack, linfo->countBlack[0]);
	fprintf(fpt, "WhiteEdge: %15.6f %d %d\n", linfo->averageWhiteEdge, linfo->sumWhiteEdge, linfo->countWhiteEdge[0]);
	fprintf(fpt, "BlackEdge: %15.6f %d %d\n", linfo->averageBlackEdge, linfo->sumBlackEdge, linfo->countBlackEdge[0]);
	fprintf(fpt, "\n");
	fprintf(fpt, "DATA: \n");
	fprintf(fpt, "%6s %10s %10s %10s %10s %15s %15s %15s %15s %15s %15s %15s %15s %15s\n", 
		"length", "White", "Black", "WhiteEdge", "BlackEdge",
		          "rateWhite", "rateBlack", "rateWhiteEdge", "rateBlackEdge",
				  "aveWhite", "aveBlack", "aveWhiteEdge", "aveBlackEdge", "realLength");
	for(i=1; i<=linfo->N; i++) {
		fprintf(fpt, "%6d %10d %10d %10d %10d %15.6f %15.6f %15.6f %15.6f %15.6f %15.6f %15.6f %15.6f %15.6f\n", 
			i,
			linfo->countWhite[i], 
			linfo->countBlack[i],
			linfo->countWhiteEdge[i], 
			linfo->countBlackEdge[i],
			(linfo->countWhite[0]?linfo->countWhite[i]/(double)linfo->countWhite[0]*100:0), 
			(linfo->countBlack[0]?linfo->countBlack[i]/(double)linfo->countBlack[0]*100:0),
			(linfo->countWhiteEdge[0]?linfo->countWhiteEdge[i]/(double)linfo->countWhiteEdge[0]*100:0), 
			(linfo->countBlackEdge[0]?linfo->countBlackEdge[i]/(double)linfo->countBlackEdge[0]*100:0),
			(linfo->sumWhite?linfo->countWhite[i]/(double)linfo->sumWhite*100:0), 
			(linfo->sumBlack?linfo->countBlack[i]/(double)linfo->sumBlack*100:0),
			(linfo->sumWhiteEdge?linfo->countWhiteEdge[i]/(double)linfo->sumWhiteEdge*100:0), 
			(linfo->sumBlackEdge?linfo->countBlackEdge[i]/(double)linfo->sumBlackEdge*100:0),
			i*linfo->Length.z);
	}
}


void
lmrcImageMeanFreePathCalcAlongZaxis(mrcImage* in, lmrcImageMeanFreePathCalcInfo* linfo, lmrcImageMeanFreePathCalcMode mode)
{
	int* countBlack;
	int* countWhite; 
	int* countBlackEdge; 
	int* countWhiteEdge; 
	int x, y, z;
	int srcx, srcy;
	//int minx, miny;
	//int maxx, maxy;
	int flagEdge;
	double data;
	double data0;
	double prevData;
	int lengthWhite;
	int lengthBlack;
	int i;
	int k;
	
	countBlack     = memoryAllocate(sizeof(int)*(in->HeaderN.z+1), "in lmrcImageMeanFreePathCalc");
	countWhite     = memoryAllocate(sizeof(int)*(in->HeaderN.z+1), "in lmrcImageMeanFreePathCalc");
	countBlackEdge = memoryAllocate(sizeof(int)*(in->HeaderN.z+1), "in lmrcImageMeanFreePathCalc");
	countWhiteEdge = memoryAllocate(sizeof(int)*(in->HeaderN.z+1), "in lmrcImageMeanFreePathCalc");

	linfo->white.Header = in->Header;	
	linfo->black.Header = in->Header;	
	linfo->whiteEdge.Header = in->Header;	
	linfo->blackEdge.Header = in->Header;	
	mrcInit(&linfo->white, NULL);
	mrcInit(&linfo->black, NULL);
	mrcInit(&linfo->whiteEdge, NULL);
	mrcInit(&linfo->blackEdge, NULL);

	linfo->countBlack = countBlack;
	linfo->countWhite = countWhite;
	linfo->countBlackEdge = countBlackEdge;
	linfo->countWhiteEdge = countWhiteEdge;
	linfo->N = in->HeaderN.z;

	for(y=0; y<in->HeaderN.y; y++) {
		linfo->min.y = floor(y - linfo->LineThickness.y/2.0 + 0.5); 
		linfo->max.y = floor(y + linfo->LineThickness.y/2.0 + 0.5); 
		linfo->cur.y = y;
	for(x=0; x<in->HeaderN.x; x++) {
		linfo->min.x = floor(x - linfo->LineThickness.x/2.0 + 0.5); 
		linfo->max.x = floor(x + linfo->LineThickness.x/2.0 + 0.5); 
		linfo->cur.x = x;

		flagEdge = 1;
		lengthWhite = 0;
		lengthBlack = 0;

		// Start
		linfo->cur.z = 0;
		data = lmrcImageMeanFreePathCalcBWEvaluation(in, linfo, mode);

		if(0<data) {
			countWhiteEdge[0]++;		
			lengthWhite++;
		} else {
			countBlackEdge[0]++;
			lengthBlack++;
		}
		prevData = data;

		// Intermediate 
		for(z=1; z<in->HeaderN.z-1; z++) {
			linfo->cur.z = z;
			data = lmrcImageMeanFreePathCalcBWEvaluation(in, linfo, mode);

			if(0<prevData) {
				if(0<data) { // White continue
					if(flagEdge) {
						countWhiteEdge[0]++;
					} else {
						countWhite[0]++;
					}
					lengthWhite++;
				} else {     // White -> Black
					if(flagEdge) {
						countWhiteEdge[lengthWhite]++;	
						for(k=z-1; z-lengthWhite<=k; k--) {
							mrcPixelDataSet(&linfo->whiteEdge, x, y, k, lengthWhite, mrcPixelRePart);
						}
						flagEdge=0;
					} else {
						countWhite[lengthWhite]++;
						for(k=z-1; z-lengthWhite<=k; k--) {
							mrcPixelDataSet(&linfo->white, x, y, k, lengthWhite, mrcPixelRePart);
						}
					}
					countBlack[0]++;
					lengthWhite=0;			
					lengthBlack++;
				}
			} else {
				if(0<data) { // Black -> White 
					if(flagEdge) {
						countBlackEdge[lengthBlack]++;	
						for(k=z-1; z-lengthBlack<=k; k--) {
							mrcPixelDataSet(&linfo->blackEdge, x, y, k, lengthBlack, mrcPixelRePart);
						}
						flagEdge=0;
					} else {
						countBlack[lengthBlack]++;
						for(k=z-1; z-lengthBlack<=k; k--) {
							mrcPixelDataSet(&linfo->black, x, y, k, lengthBlack, mrcPixelRePart);
						}
					}
					countWhite[0]++;
					lengthBlack=0;
					lengthWhite++;
				} else {    // Black continue
					if(flagEdge) {
						countBlackEdge[0]++;	
					} else {
						countBlack[0]++;	
					}
					lengthBlack++;			
				}
			}
			prevData = data;
		}

		// End
		z = in->HeaderN.z-1;
		linfo->cur.z = z;
		data = lmrcImageMeanFreePathCalcBWEvaluation(in, linfo, mode);

		if(0<prevData) {
			if(0<data) { // White continue
				lengthWhite++;
				for(k=z; z-lengthWhite<k; k--) {
					mrcPixelDataSet(&linfo->whiteEdge, x, y, k, lengthWhite, mrcPixelRePart);
				}
				countWhiteEdge[lengthWhite]++;	
				countWhiteEdge[0]++;
			} else {     // White -> Black
				if(flagEdge) {
					countWhiteEdge[lengthWhite]++;
					for(k=z-1; z-lengthWhite<=k; k--) {
						mrcPixelDataSet(&linfo->whiteEdge, x, y, k, lengthWhite, mrcPixelRePart);
					}
				} else {
					countWhite[lengthWhite]++;
					for(k=z-1; z-lengthWhite<=k; k--) {
						mrcPixelDataSet(&linfo->white, x, y, k, lengthWhite, mrcPixelRePart);
					}
				}
				mrcPixelDataSet(&linfo->blackEdge, x, y, z, 1, mrcPixelRePart);
				countBlackEdge[1]++;	
				countBlackEdge[0]++;
			}
		} else {
			if(0<data) { // Black -> White 
				if(flagEdge) {
					countBlackEdge[lengthBlack]++;	
					for(k=z-1; z-lengthBlack<=k; k--) {
						mrcPixelDataSet(&linfo->blackEdge, x, y, k, lengthBlack, mrcPixelRePart);
					}
				} else {
					countBlack[lengthBlack]++;
					for(k=z-1; z-lengthBlack<=k; k--) {
						mrcPixelDataSet(&linfo->black, x, y, k, lengthBlack, mrcPixelRePart);
					}
				}
				mrcPixelDataSet(&linfo->whiteEdge, x, y, z, 1, mrcPixelRePart);
				countWhiteEdge[1]++;	
				countWhiteEdge[0]++;
			} else {    // Black continue
				lengthBlack++;
				countBlackEdge[lengthBlack]++;	
				for(k=z; z-lengthBlack<k; k--) {
					mrcPixelDataSet(&linfo->blackEdge, x, y, k, lengthBlack, mrcPixelRePart);
				}
				countBlackEdge[0]++;	
			}
		}
	}
	}

	linfo->averageWhite = linfo->averageBlack = 0; 
	linfo->averageWhiteEdge = linfo->averageBlackEdge = 0; 
	linfo->sumWhite = linfo->sumBlack = 0; 
	linfo->sumWhiteEdge = linfo->sumBlackEdge = 0; 
	for(i=1; i<=linfo->N; i++) {
		linfo->averageWhite	    += i*i*linfo->countWhite[i];		
		linfo->averageBlack     += i*i*linfo->countBlack[i];		
		linfo->averageWhiteEdge	+= i*i*linfo->countWhiteEdge[i];			
		linfo->averageBlackEdge += i*i*linfo->countBlackEdge[i];			
		linfo->sumWhite	    += linfo->countWhite[i];		
		linfo->sumBlack     += linfo->countBlack[i];		
		linfo->sumWhiteEdge	+= linfo->countWhiteEdge[i];			
		linfo->sumBlackEdge += linfo->countBlackEdge[i];			
	}

	if(0<linfo->countWhite[0]) {
		linfo->averageWhite	    /= linfo->countWhite[0];		
	} else { 
		linfo->averageWhite	    = 0;
	}
	if(0<linfo->countBlack[0]) {
		linfo->averageBlack	    /= linfo->countBlack[0];		
	} else { 
		linfo->averageBlack	    = 0;
	}
	if(0<linfo->countWhiteEdge[0]) {
		linfo->averageWhiteEdge	    /= linfo->countWhiteEdge[0];		
	} else { 
		linfo->averageWhiteEdge	    = 0;
	}
	if(0<linfo->countBlackEdge[0]) {
		linfo->averageBlackEdge	    /= linfo->countBlackEdge[0];		
	} else { 
		linfo->averageBlackEdge	    = 0;
	}

	if(MeanFreePathCalcModeRealLength&mode) {
		DEBUGPRINT1("mode: %d in Calc: RealLength Start\n", mode);
		lmrcImageMultiplying2(&linfo->white, linfo->white.HeaderLength.z, linfo->white.HeaderMode);
		lmrcImageMultiplying2(&linfo->whiteEdge, linfo->whiteEdge.HeaderLength.z, linfo->whiteEdge.HeaderMode);
		lmrcImageMultiplying2(&linfo->black, linfo->black.HeaderLength.z, linfo->black.HeaderMode);
		lmrcImageMultiplying2(&linfo->blackEdge, linfo->blackEdge.HeaderLength.z, linfo->blackEdge.HeaderMode);
		DEBUGPRINT1("mode: %d in Calc: RealLength End\n", mode);
	}
}


double 
lmrcImageMeanFreePathCalcBWEvaluation(mrcImage* in, lmrcImageMeanFreePathCalcInfo* linfo, lmrcImageMeanFreePathCalcMode  mode)
{
	double data, data0;
	double rx, ry, r0x, r0y;
	mrcImageParaTypeReal srcx, srcy;
	int max;

	switch(linfo->LineShape) {
		case lmrcImageMeanFreePathCalcLineShapeRectangle: {
			switch(linfo->LineMode) {
				case lmrcImageMeanFreePathCalcLineModeAnd: {
					data = 1;
					for(srcy=linfo->min.y; srcy<linfo->max.y; srcy++) { 
					for(srcx=linfo->min.x; srcx<linfo->max.x; srcx++) { 
						mrcPixelDataGet(in, srcx, srcy, linfo->cur.z, &data0, mrcPixelRePart, mrcPixelHowNearest);
						data *= data0;
					}
					}
					if(0<data) data = 1;
					break;
				}
				case lmrcImageMeanFreePathCalcLineModeOr: {
					data = 0;
					for(srcy=linfo->min.y; srcy<linfo->max.y; srcy++) { 
					for(srcx=linfo->min.x; srcx<linfo->max.x; srcx++) { 
						mrcPixelDataGet(in, srcx, srcy, linfo->cur.z, &data0, mrcPixelRePart, mrcPixelHowNearest);
						data += data0;
						if(0<data) break;
					}
						if(0<data) break;
					}
					if(0<data) data = 1;
					break;
				}
				case lmrcImageMeanFreePathCalcLineModeGE: {
					data = 0;
					max  = 0;
					for(srcy=linfo->min.y; srcy<linfo->max.y; srcy++) { 
					for(srcx=linfo->min.x; srcx<linfo->max.x; srcx++) { 
						mrcPixelDataGet(in, srcx, srcy, linfo->cur.z, &data0, mrcPixelRePart, mrcPixelHowNearest);
						data += data0;
						max++;
					}
					}
					if(max/2<=data) {
						data = 1;
					}
					break;
				}
				default: {
					fprintf(stderr, "Not supported Line Mode: %d\n", linfo->LineMode);
					exit(EXIT_FAILURE);
					break;
				}
			}
			break;
		}
		case lmrcImageMeanFreePathCalcLineShapeDisk: {
			r0x = (linfo->max.x - linfo->cur.x);
			r0y = (linfo->max.y - linfo->cur.y);
			switch(linfo->LineMode) {
				case lmrcImageMeanFreePathCalcLineModeAnd: {
					data = 1;
					for(srcy=linfo->min.y; srcy<=linfo->max.y; srcy++) { 
						ry = SQR((srcy - linfo->cur.y)/r0y); 
					for(srcx=linfo->min.x; srcx<linfo->max.x; srcx++) { 
						rx = SQR((srcx - linfo->cur.x)/r0x); 

						if(rx + ry <= 1) {
							mrcPixelDataGet(in, srcx, srcy, linfo->cur.z, &data0, mrcPixelRePart, mrcPixelHowNearest);
							data *= data0;
						}
					}
					}
					if(0<data) data = 1;
					break;
				}
				case lmrcImageMeanFreePathCalcLineModeOr: {
					data = 0;
					for(srcy=linfo->min.y; srcy<linfo->max.y; srcy++) { 
						ry = SQR((srcy - linfo->cur.y)/r0y); 
					for(srcx=linfo->min.x; srcx<linfo->max.x; srcx++) { 
						rx = SQR((srcx - linfo->cur.x)/r0x); 

						if(rx + ry <= 1) {
							mrcPixelDataGet(in, srcx, srcy, linfo->cur.z, &data0, mrcPixelRePart, mrcPixelHowNearest);
							data += data0;
							if(0<data) break;
						}
					}
						if(0<data) break;
					}
					if(0<data) data = 1;
					break;
				}
				case lmrcImageMeanFreePathCalcLineModeGE: {
					data = 0;
					max  = 0;
					for(srcy=linfo->min.y; srcy<linfo->max.y; srcy++) { 
						ry = SQR((srcy - linfo->cur.y)/r0y); 
					for(srcx=linfo->min.x; srcx<linfo->max.x; srcx++) { 
						rx = SQR((srcx - linfo->cur.x)/r0x); 

						if(rx + ry <= 1) {
							mrcPixelDataGet(in, srcx, srcy, linfo->cur.z, &data0, mrcPixelRePart, mrcPixelHowNearest);
							data += data0;
							max++;
						}
					}
					}
					if(max/2<=data) {
						data = 1;
					}
					break;
				}
				default: {
					fprintf(stderr, "Not supported Line Mode: %d\n", linfo->LineMode);
					exit(EXIT_FAILURE);
					break;
				}
			}
			break;
		}
		default: {
			fprintf(stderr, "Not supported Line Shape: %d\n", linfo->LineShape);
			exit(EXIT_FAILURE);
			break;
		}
	}
	return data;
}

