1 /*
2 commonLinesReg -- Peter Lauren, 2020-08-21
3
4 SYNOPSIS
5
6 commonLinesReg -i <Input Filename> -p <ac|as|cs|acs> -o <orientation> [-f [-l <min. length>]]
7
8 DESCRIPTION
9
10 Uses common lines algorithm to linearly align two volumes. It does this by finding where orthogonal volumes intersect
11 in Fourier space. If they are rotationally aligned, they should intersect on the appropriate orthogonal cardinal axes.
12 The slope of the phase shift, along the lines of intersection, reflect the translational misalignment. The translations,
13 and rotations, are each determined by solving a linear system. These solutions are completely independent of each other.
14
15 -o Output of "3dinfo -orient mydset" in the linux command line. This gives the a, y,z z-axes in terms of the axial, coronal
16 and sagmital direction.
17
18 -f Examine effects of frequency range. The search is over a range of spectral lengths and min. frequencies. "min. length"
19 is the starting length.
20
21 NOTES
22
23 1/ The current implementation aligns volumes made from unweaving the interlaced sections of a brick.
24 2/ Two to three projection axes must be chosen:
25 - ac: Axial and Coronal
26 - as: Axial and Sagital
27 - cs: Coronal and Sagital
28 - acs: Axial, Coronal and Sagital
29
30 */
31
32 #include <stdio.h>
33 #include <dirent.h>
34 #include <malloc.h>
35 #include <string.h>
36 #include <unistd.h>
37 #include <sys/resource.h>
38 #include "mrilib.h"
39 #include "Generic.h"
40 // #include "VolumeRegistration.h"
41
42 #define UNIX 1
43
44 int Cleanup(char *inputFileName, COMPLEX ***TwoDFts, THD_3dim_dataset *din);
45 int open_input_dset(THD_3dim_dataset ** din, char * fname);
46 int unweave_sections(THD_3dim_dataset *din, THD_3dim_dataset **dodd, THD_3dim_dataset **deven);
47 int makeProjection(THD_3dim_dataset *din, THD_3dim_dataset **dout, char projCode, char *orientation);
48 int float2DImage(THD_3dim_dataset *dset);
49 int getLargestDimension(THD_3dim_dataset *din);
50 int zeroPadProjection(THD_3dim_dataset *din, THD_3dim_dataset **dout, int largestDimension, int factor);
51 int doesFileExist(char * searchPath, char * prefix,char *appendage , char * outputFileName);
52 int FFT2D(COMPLEX **c,int nx,int ny,int dir);
53 int FFT(int dir,int m,double *x,double *y);
54 int Powerof2(int n,int *m,int *twopm);
55 int get2DFourierTransform(THD_3dim_dataset *din, COMPLEX ***TwoDFts);
56 int outputFourierSpectrum(COMPLEX **TwoDFts, THD_3dim_dataset *din);
57 int outputRealAndImaginaryPlanes(COMPLEX **TwoDFts, THD_3dim_dataset *din);
58 COMPLEX **CheckeredInversion(COMPLEX **input, unsigned int dimension);
59 ERROR_NUMBER OutputFourierIntersectionMap(COMPLEX ***TwoDFts, int lNumberOfImages, int dimension,
60 IntRange irFrequencyRange, char *csSearchPath, char *csRootName, float *fpCoplanarAngularShift);
61 ERROR_NUMBER MakeRadialPhaseSampleArray(ComplexPlane cpFourierTransform, float fIncrementInDegrees,
62 IntRange lrFrequencyRange, float ***fpppRadialSamples, long *lpNumRadialSamples);
63 ERROR_NUMBER RotScaleComplexPlaneM(ComplexPlane *cppComplexPlane, Matrix mTransMatrix, BOOL bClip, BYTE bInterpOrder);
64 ERROR_NUMBER RotScaleComplexPlane(ComplexPlane *cppComplexPlane, FloatVector fvTranslation,
65 float fClockwiseRotationInDegrees, float fScaleX, float fScaleY, BOOL bClip, BYTE bInterpOrder);
66 ERROR_NUMBER MakeCardinalRadialPhaseSamples(FloatPlane fpReal, FloatPlane fpImaginary,
67 short *spCardinalIndices, IntRange lrFrequencyRange, float ***fpppRadialSamples);
68 ERROR_NUMBER ComparePhaseSamples(float *fpTargetSample, float *fpRefSample,
69 IntRange lrFrequencyRange, float *fpComparisonMetric);
70 float ExpectedAngularShift(int iCurrentAxis);
71 float GetAngularShift(FloatVector fvCoords);
72 void DetermineUnclippedOutputDimensions(long lInputColumns, long lInputRows,
73 float fClockwiseRotationInDegrees, FloatVector fvTranslation, float fScaleX, float fScaleY);
74 void MatVectMult(Matrix mMatrix, float *fpInVector, float *fpOutVector);
75 ERROR_NUMBER MakeRadialPhaseVectorSample(COMPLEX *cpComplexRadialSample, int iLength, float **fppRadialSample);
76 ERROR_NUMBER MakePhaseShiftVector(float *fpTargetSample, float *fpRefSample, IntRange lrFrequencyRange,
77 long *lpVectorLength, float **fppPhaseShiftVector);
78 void GetSlope(float *fVector, IntRange lrFrequencyRange, float *fpSlope, float *fpIntercept);
79 float GetMSE(float *fVector, float fSlope, float fPhaseShiftIntercept, IntRange lrFrequencyRange);
80 ERROR_NUMBER MedianFilterFloatVector(float *fpVector, long lLength);
81 float Float3x1Median(float *fpInput);
82 void FreeFloatPlane(FloatPlane Float_Plane);
83 int Round(float arg);
84 ERROR_NUMBER ErrorOpeningFile(char *csFileName);
85 float FloatPythag(float fArg1, float fArg2);
86 double Pythag(double dArg1, double dArg2);
87 double Interp_BL(float **image, float x_coord, float y_coord);
88 double Interp_BQ(float **image, float x_coord, float y_coord);
89 double Interp_BC(float **image, float x_coord, float y_coord);
90 long DRound(double arg);
91 ERROR_NUMBER Complex2Float(ComplexPlane cpInputPlane, FloatPlane *fppRealOut, FloatPlane *fppImagOut);
92 FloatPlane MakeFloatPlane(long iRows, long iColumns, float **Data);
93 char * RootName(char * csFullPathName);
94 ERROR_NUMBER GetDirectory(char *csInputFileName, char *csDirectory);
95 ERROR_NUMBER shortToFloat(THD_3dim_dataset **din);
96 ERROR_NUMBER analyzeFrequencyRange(COMPLEX ***TwoDFts, int numberOfImages,
97 int paddedDimension, char *searchPath, char *prefix, int minLength, int startingTargetIndex);
98 ERROR_NUMBER getImsePeakAndMean(COMPLEX ***TwoDFts, int dimension, int refIndex, int targetIndex,
99 IntRange irFrequencyRange, float *maxIMSE, float *meanIMSE);
100
101
102 double dCommonLinesAngleErr;
103 float fCommonLinesPeak, fCommonLinesPeakDiff;
104 LongVector lvOutputDimensions;
105
106
main(int argc,char * argv[])107 int main( int argc, char *argv[] ) {
108 ERROR_NUMBER enErrorNumber;
109 IntRange irFrequencyRange;
110 float fCoplanarAngularShift;
111 THD_3dim_dataset * din = NULL, *dodd, *deven;
112 THD_3dim_dataset* projections[6], *paddedProjections[6];
113 int i, j, largestDimension, paddingFactor=4, minLength=8;
114 char *inputFileName=NULL, projectionString[3]={'a','s','\0'};
115 char orientation[8];
116 COMPLEX** TwoDFts[6]={NULL, NULL, NULL};
117 Boolean frequencyRangeEffects = false;
118 int startingTargetIndex=4;
119
120 sprintf(orientation, "ASL");
121
122 paddingFactor=2; // DEBUG
123
124 for (i=0; i<argc; ++i) if (argv[i][0]=='-'){
125 switch(argv[i][1]){
126 case 'f':
127 frequencyRangeEffects = true;
128 break;
129
130 case 'i':
131 if (!(inputFileName=(char*)malloc(strlen(argv[++i])+8)))
132 return Cleanup(inputFileName, TwoDFts, din);
133 sprintf(inputFileName,"%s",argv[i]);
134 break;
135
136 case 'l':
137 minLength=atoi(argv[++i]);
138 break;
139
140 case 'o':
141 sprintf(orientation, "%s", argv[++i]);
142 break;
143
144 case 'p':
145 sprintf(projectionString,"%s",argv[++i]);
146 break;
147
148 case 's':
149 startingTargetIndex=atoi(argv[++i]);
150 break;
151 }
152 }
153
154 if( open_input_dset(&din, inputFileName) )
155 return Cleanup(inputFileName, TwoDFts, din);
156
157 // Split volume into odd and even slices
158 if (!unweave_sections(din, &dodd, &deven))
159 return Cleanup(inputFileName, TwoDFts, din);
160
161 // Make odd projections
162 makeProjection(dodd, &(projections[0]), 's', orientation);
163 makeProjection(dodd, &(projections[1]), 'c', orientation);
164 makeProjection(dodd, &(projections[2]), 'a', orientation);
165
166 // Make even projections
167 makeProjection(deven, &(projections[3]), 's', orientation);
168 makeProjection(deven, &(projections[4]), 'c', orientation);
169 makeProjection(deven, &(projections[5]), 'a', orientation);
170
171 // Get largest dimension from data set
172 largestDimension=getLargestDimension(dodd);
173
174 // Free up intput dataset
175 DSET_delete(din);
176
177 // Float projections
178 // Currently skipped because of head voxels on the edge of the projection image. Apart from that,
179 // the edge defaults to zero. Floating creates edge effects rather than suppressing them.
180 // float2DImage(doddSagProj);
181
182 // Zero pad projections
183 for (i=0; i<6; ++i){
184 if (!(zeroPadProjection(projections[i], &(paddedProjections[i]), largestDimension, paddingFactor))){
185 Cleanup(inputFileName, TwoDFts, projections[2]);
186 for (j=0; j<i; ++j) DSET_delete(paddedProjections[j]);
187 for (; j<6; ++j) DSET_delete(projections[j]);
188 }
189 DSET_delete(projections[i]);
190 }
191
192 // Make Fourier transforms of projections
193 for (i=0; i<6; ++i){
194 if (!(get2DFourierTransform(paddedProjections[i], &(TwoDFts[i])))){
195 Cleanup(inputFileName, TwoDFts, paddedProjections[0]);
196 for (i=0; i<6; ++i) DSET_delete(paddedProjections[i]);
197 }
198 }
199
200 // Output Fourier spectra, of projections.
201 #if 1
202 for (i=0; i<6; ++i){
203 if (!(outputFourierSpectrum(TwoDFts[i], paddedProjections[i]))){
204 Cleanup(inputFileName, TwoDFts, paddedProjections[0]);
205 for (i=0; i<6; ++i) DSET_delete(paddedProjections[i]);
206 }
207 }
208 #endif
209
210 // DEBUG: Output real and imaginary plane of first FT.
211 // outputRealAndImaginaryPlanes(TwoDFts[0], paddedProjections[0]);
212
213 // Make CSV file of radial phase shift linearity (RPSL) between pairs of FTs
214 int paddedDimension = DSET_NY(paddedProjections[0]);
215 int lNumberOfImages = 6;
216 char *prefix=DSET_PREFIX(din);
217 char *searchPath=(char *)malloc(strlen(inputFileName)*sizeof(char));
218 if (enErrorNumber=GetDirectory(inputFileName, searchPath)){
219 Cleanup(inputFileName, TwoDFts, paddedProjections[0]);
220 for (i=0; i<6; ++i) DSET_delete(paddedProjections[i]);
221 free(searchPath);
222 return enErrorNumber;
223 }
224 if (frequencyRangeEffects){
225 if ((enErrorNumber=analyzeFrequencyRange(TwoDFts, lNumberOfImages, paddedDimension,
226 searchPath, prefix, minLength, startingTargetIndex))!=ERROR_NONE){
227 Cleanup(inputFileName, TwoDFts, paddedProjections[0]);
228 for (i=0; i<6; ++i) DSET_delete(paddedProjections[i]);
229 free(searchPath);
230 return enErrorNumber;
231 }
232 } else {
233 irFrequencyRange.iMin = 15;
234 irFrequencyRange.iMax = paddedDimension/2-32;
235 if ((enErrorNumber=OutputFourierIntersectionMap(TwoDFts, lNumberOfImages, paddedDimension,
236 irFrequencyRange, searchPath, prefix, &fCoplanarAngularShift))!=ERROR_NONE){
237 Cleanup(inputFileName, TwoDFts, paddedProjections[0]);
238 for (i=0; i<6; ++i) DSET_delete(paddedProjections[i]);
239 free(searchPath);
240 return enErrorNumber;
241 }
242 }
243 free(searchPath);
244
245 Cleanup(inputFileName, TwoDFts, paddedProjections[0]);
246 for (i=0; i<6; ++i) DSET_delete(paddedProjections[i]);
247
248 return 1;
249 }
250
MakeRadialPhaseSampleArray(ComplexPlane cpFourierTransform,float fIncrementInDegrees,IntRange lrFrequencyRange,float *** fpppRadialSamples,long * lpNumRadialSamples)251 ERROR_NUMBER MakeRadialPhaseSampleArray(ComplexPlane cpFourierTransform, float fIncrementInDegrees,
252 IntRange lrFrequencyRange, float ***fpppRadialSamples, long *lpNumRadialSamples)
253 {
254 ERROR_NUMBER enErrorNumber;
255 long iSampleIndex, iQuadIndex;
256 float fAngleInDegrees=0;
257 ComplexPlane cpRotatedFT;
258 FloatPlane fpReal, fpImag;
259 FloatVector fvTranslationVector={0,0};
260 short saCardinalIndices[4];
261
262 /* TODO: Add later
263 if (bSincFunctionInterp) return MakeRadialPhaseSampleArrayWithSincInterp(fpRealInput, fIncrementInDegrees,
264 lrFrequencyRange, fpppRadialSamples, lpNumRadialSamples, sFFTPaddingFactor);
265 */
266
267 // Make rotation buffer
268 cpRotatedFT=MakeComplexPlane(cpFourierTransform.iRows, cpFourierTransform.iColumns, NULL);
269 if (!(cpRotatedFT.Data)) return ERROR_MEMORY_ALLOCATION;
270
271 // Allocate pointer to array
272 if (!(*fpppRadialSamples=(float **)malloc(4*((*lpNumRadialSamples)+4)*sizeof(float *))))
273 {
274 FreeComplexPlane(cpRotatedFT);
275 return ERROR_MEMORY_ALLOCATION;
276 }
277
278 // Initlialize cardinal indices
279 saCardinalIndices[0]=0;
280 saCardinalIndices[1]=(short)(*lpNumRadialSamples);
281 saCardinalIndices[2]=(short)((*lpNumRadialSamples)*2);
282 saCardinalIndices[3]=saCardinalIndices[1]+saCardinalIndices[2];
283
284 // Make samples
285 for (iSampleIndex=0; iSampleIndex<(*lpNumRadialSamples); ++iSampleIndex)
286 {
287 // Copy input image into rotation buffer
288 CopyComplexPlane(cpFourierTransform, cpRotatedFT);
289
290 // Apply rotation to spatial plane image
291 if (fAngleInDegrees!=0 && (enErrorNumber=RotScaleComplexPlane(&cpRotatedFT,
292 fvTranslationVector, -fAngleInDegrees, 1.0f, 1.0f, TRUE, BIQUADRATIC_INTERPOLATION))!=ERROR_NONE)
293 {
294 FreeComplexPlane(cpRotatedFT);
295 for (--iSampleIndex; iSampleIndex>=0; --iSampleIndex) free((*fpppRadialSamples)[iSampleIndex]);
296 free(*fpppRadialSamples);
297 (*fpppRadialSamples)=NULL;
298 return enErrorNumber;
299 }
300
301 // Make real and imaginary planes from Fourier transform
302 if ((enErrorNumber=Complex2Float(cpRotatedFT, &fpReal, &fpImag))!=ERROR_NONE)
303 {
304 FreeComplexPlane(cpRotatedFT);
305 for (--iSampleIndex; iSampleIndex>=0; --iSampleIndex) free((*fpppRadialSamples)[iSampleIndex]);
306 free(*fpppRadialSamples);
307 (*fpppRadialSamples)=NULL;
308 return enErrorNumber;
309 }
310
311 // Make radial samples from cardinal angles
312 if ((enErrorNumber=MakeCardinalRadialPhaseSamples(fpReal, fpImag, &(saCardinalIndices[0]),
313 lrFrequencyRange, fpppRadialSamples))!=ERROR_NONE)
314 {
315 FreeComplexPlane(cpRotatedFT);
316 FreeFloatPlane(fpReal);
317 FreeFloatPlane(fpImag);
318 for (--iSampleIndex; iSampleIndex>=0; --iSampleIndex) free((*fpppRadialSamples)[iSampleIndex]);
319 free ((*fpppRadialSamples));
320 (*fpppRadialSamples)=NULL;
321 return enErrorNumber;
322 }
323
324
325 // Increment angle and cardinal indices
326 fAngleInDegrees+=fIncrementInDegrees;
327 for (iQuadIndex=0; iQuadIndex<4; ++iQuadIndex) ++saCardinalIndices[iQuadIndex];
328
329 // Cleanup
330 FreeFloatPlane(fpReal);
331 FreeFloatPlane(fpImag);
332 }
333
334 // Cleanup
335 FreeComplexPlane(cpRotatedFT);
336
337 return ERROR_NONE;
338 }
339
analyzeFrequencyRange(COMPLEX *** TwoDFts,int numberOfImages,int paddedDimension,char * searchPath,char * prefix,int minLength,int startingTargetIndex)340 ERROR_NUMBER analyzeFrequencyRange(COMPLEX ***TwoDFts, int numberOfImages,
341 int paddedDimension, char *searchPath, char *prefix, int minLength, int startingTargetIndex){
342
343 ERROR_NUMBER enErrorNumber;
344 char csAnalysisFileName[512];
345 IntRange irFrequencyRange;
346 FILE *outputFile;
347 int maxMin=paddedDimension/4, maxMax=paddedDimension/2-1;
348 int minRangeLength=minLength;
349 // int minRangeLength=250; // DEBUG
350 float maxIMSE, meanIMSE; // Max and mean inverse mean squared error
351 struct rusage r_usage; // Track memory usage
352
353 for (int lPlaneIndex=startingTargetIndex; lPlaneIndex<numberOfImages; ++lPlaneIndex){ // Avoid coplanar and self reference
354 fprintf(stdout, "Processing image %d\n", lPlaneIndex);
355
356 // Open output file
357 sprintf(csAnalysisFileName,"%s/%s_LOI_freqRange%d.csv", searchPath, prefix, lPlaneIndex);
358 if (!(outputFile=fopen(csAnalysisFileName, "w")))
359 return ErrorOpeningFile(csAnalysisFileName);
360
361 // Make header
362 fprintf(outputFile, "Range Length\tFMin\n");
363 for (int i=0; i<=maxMin; ++i) fprintf(outputFile, "\t%d", i);
364 fprintf(outputFile, "\n");
365
366 // Process for each frequency range length
367 for (int rangeLength=minRangeLength; rangeLength<maxMax; ++rangeLength){
368 getrusage(RUSAGE_SELF,&r_usage);
369 fprintf(stderr, "Processing length %d of %d. Memory usage: %ld kilobytes\n",
370 rangeLength, maxMax, r_usage.ru_maxrss);
371 fprintf(outputFile, "%d", rangeLength);
372
373 for (int fMin=0; fMin<=maxMin; ++fMin){
374 fprintf(stderr, "fMin=%d of %d\r", fMin, maxMin);
375 if (fMin+rangeLength>maxMax) fprintf(outputFile, "\t0");
376 else {
377 irFrequencyRange.iMin = fMin;
378 irFrequencyRange.iMax = fMin+rangeLength;
379
380 if ((enErrorNumber=getImsePeakAndMean(TwoDFts, paddedDimension, 0, lPlaneIndex,
381 irFrequencyRange, &maxIMSE, &meanIMSE))!=ERROR_NONE){
382 fclose(outputFile);
383 return enErrorNumber;
384 }
385
386 fprintf(outputFile, "\t%f", maxIMSE/meanIMSE);
387 }
388 }
389 fprintf(outputFile, "\n");
390 fprintf(stderr, "\n");
391 }
392
393 fclose(outputFile);
394 }
395
396 return ERROR_NONE;
397 }
398
getImsePeakAndMean(COMPLEX *** TwoDFts,int dimension,int refIndex,int targetIndex,IntRange irFrequencyRange,float * maxIMSE,float * meanIMSE)399 ERROR_NUMBER getImsePeakAndMean(COMPLEX ***TwoDFts, int dimension, int refIndex, int targetIndex, IntRange irFrequencyRange,
400 float *maxIMSE, float *meanIMSE){
401
402 ERROR_NUMBER enErrorNumber;
403 ComplexPlane cpComplexPlane;
404 float fIncrementInDegrees=2.0f*(float)(asin(0.5/irFrequencyRange.iMax)/DEGREES2RADIANS);
405 float **fppRefSamples, **fpTargetSamples, fComparisonMetric;
406 long lPlaneIndex, lNumTargetSamples, lNumRefSamples, lTargetIndex, lRefIndex;
407
408 // Make reference samples
409 lNumRefSamples=Round(180.0f/fIncrementInDegrees)-1;
410 cpComplexPlane.iColumns=cpComplexPlane.iRows=dimension;
411 cpComplexPlane.Data=TwoDFts[refIndex];
412 if ((enErrorNumber=MakeRadialPhaseSampleArray(cpComplexPlane, fIncrementInDegrees,
413 irFrequencyRange, &fppRefSamples, &lNumRefSamples))!=ERROR_NONE){
414 return ERROR_MEMORY_ALLOCATION;
415 }
416
417 // Make target array
418 lNumTargetSamples=Round(360.0f/fIncrementInDegrees);
419 cpComplexPlane.iColumns=cpComplexPlane.iRows=dimension;
420 cpComplexPlane.Data=TwoDFts[targetIndex];
421 if ((enErrorNumber=MakeRadialPhaseSampleArray(cpComplexPlane, fIncrementInDegrees,
422 irFrequencyRange, &fpTargetSamples, &lNumTargetSamples))!=ERROR_NONE)
423 {
424 for (lRefIndex=0; lRefIndex<lNumRefSamples; ++lRefIndex)
425 free(fppRefSamples[lRefIndex]);
426 free(fppRefSamples);
427 return enErrorNumber;
428 }
429
430 // Process samples from reference and target plane
431 *maxIMSE = *meanIMSE = 0;
432 for (lTargetIndex=0; lTargetIndex<=lNumTargetSamples; ++lTargetIndex)
433 {
434 // Compare target samples against each reference sample using multiplicative inverse of
435 // mean square distance (MSD)
436 for (lRefIndex=0; lRefIndex<lNumRefSamples; ++lRefIndex)
437 {
438 // Get comparison metric between the two samples
439 if ((enErrorNumber=ComparePhaseSamples(fpTargetSamples[lTargetIndex], fppRefSamples[lRefIndex],
440 irFrequencyRange, &fComparisonMetric))!=ERROR_NONE)
441 {
442 for (lRefIndex=0; lRefIndex<lNumRefSamples; ++lRefIndex) free(fppRefSamples[lRefIndex]);
443 free(fppRefSamples);
444 for (lTargetIndex=0; lTargetIndex<lNumRefSamples; ++lTargetIndex) free(fpTargetSamples[lTargetIndex]);
445 free(fpTargetSamples);
446 return enErrorNumber;
447 }
448
449 // Update mean
450 *meanIMSE += fComparisonMetric;
451
452 // Update max
453 if (fComparisonMetric>*maxIMSE) *maxIMSE = fComparisonMetric;
454 }
455 }
456
457 *meanIMSE /= lNumTargetSamples*lNumRefSamples;
458
459 // Cleanup
460 int numberToFree = lNumRefSamples*4;
461 for (lRefIndex=0; lRefIndex<numberToFree; ++lRefIndex) free(fppRefSamples[lRefIndex]);
462 free(fppRefSamples);
463 numberToFree = lNumTargetSamples*4;
464 for (lTargetIndex=0; lTargetIndex<numberToFree; ++lTargetIndex) free(fpTargetSamples[lTargetIndex]);
465 free(fpTargetSamples);
466
467 return ERROR_NONE;
468 }
469
OutputFourierIntersectionMap(COMPLEX *** TwoDFts,int lNumberOfImages,int dimension,IntRange irFrequencyRange,char * csSearchPath,char * csRootName,float * fpCoplanarAngularShift)470 ERROR_NUMBER OutputFourierIntersectionMap(COMPLEX ***TwoDFts, int lNumberOfImages, int dimension,
471 IntRange irFrequencyRange, char *csSearchPath, char *csRootName, float *fpCoplanarAngularShift)
472 {
473 ERROR_NUMBER enErrorNumber;
474 ComplexPlane cpComplexPlane;
475 char csAnalysisFileName[512];
476
477 // Increment to avoid false maxima
478 float fIncrementInDegrees=2.0f*(float)(asin(0.5/irFrequencyRange.iMax)/DEGREES2RADIANS);
479
480 float **fppRefSamples, **fpTargetSamples, fComparisonMetric, faPeakValues[10];
481 long lPlaneIndex, lNumTargetSamples, lNumRefSamples, lTargetIndex, lRefIndex;
482 short sPeakIndex, i, j;
483 FloatVector fvaPeakCoords[10];
484 FILE *fpAnalysisFile, *fpGraphFile;
485 BOOL bAnalysisFile=TRUE;
486 char *graphFileName;
487 // , *directory, *rootName;
488
489 // Try reducing increment to improve accuracy
490 fIncrementInDegrees /= 2;
491
492 // Allocate memory to graph file name
493 if (!(graphFileName=(char *)malloc(strlen(csAnalysisFileName)+32)))
494 return ERROR_MEMORY_ALLOCATION;
495
496 // Make reference samples
497 lNumRefSamples=Round(180.0f/fIncrementInDegrees)-1;
498 cpComplexPlane.iColumns=cpComplexPlane.iRows=dimension;
499 cpComplexPlane.Data=TwoDFts[0];
500 if ((enErrorNumber=MakeRadialPhaseSampleArray(cpComplexPlane, fIncrementInDegrees,
501 irFrequencyRange, &fppRefSamples, &lNumRefSamples))!=ERROR_NONE){
502 free(graphFileName);
503 return ERROR_MEMORY_ALLOCATION;
504 }
505
506 /*
507 // Get directory name
508 if ((enErrorNumber=GetDirectory(csAnalysisFileName, directory))!=ERROR_NONE){
509 free(graphFileName);
510 free(directory);
511 free(rootName);
512 return enErrorNumber;
513 }
514 */
515
516 // Process target planes
517 for (lPlaneIndex=1; lPlaneIndex<lNumberOfImages; ++lPlaneIndex)
518 if (lPlaneIndex!=3) // Avoid coplanar
519 {
520 // Open graph file
521 sprintf(graphFileName, "%s/%s_LOI_graph%ld.csv", csSearchPath, csRootName, lPlaneIndex);
522 if (!(fpGraphFile=fopen(graphFileName, "w"))){
523 free(graphFileName);
524 return ErrorOpeningFile(graphFileName);
525 }
526
527 // Make target array
528 lNumTargetSamples=Round(360.0f/fIncrementInDegrees);
529 cpComplexPlane.iColumns=cpComplexPlane.iRows=dimension;
530 cpComplexPlane.Data=TwoDFts[lPlaneIndex];
531 if ((enErrorNumber=MakeRadialPhaseSampleArray(cpComplexPlane, fIncrementInDegrees,
532 irFrequencyRange, &fpTargetSamples, &lNumTargetSamples))!=ERROR_NONE)
533 {
534 for (lRefIndex=0; lRefIndex<lNumRefSamples; ++lRefIndex)
535 free(fppRefSamples[lRefIndex]);
536 fclose(fpGraphFile);
537 free(fppRefSamples);
538 free(graphFileName);
539 return enErrorNumber;
540 }
541
542 // Write graphic file header
543 fprintf(fpGraphFile, "Ref.\\Target");
544 for (lRefIndex=0; lRefIndex<=lNumRefSamples; ++lRefIndex)
545 fprintf(fpGraphFile, ",%3.2f", fIncrementInDegrees*lRefIndex);
546 fprintf(fpGraphFile, "\n");
547
548 // Initialize peak values (Maybe should go in plane loop)
549 for (i=0; i<10; ++i) faPeakValues[i]=0.0f;
550
551 // Process samples from target plane
552 for (lTargetIndex=0; lTargetIndex<=lNumTargetSamples; ++lTargetIndex)
553 {
554 // Write graph file entry
555 fprintf(fpGraphFile, "%3.2f", fIncrementInDegrees*lTargetIndex);
556
557 // Compare target samples against each reference sample using multiplicative inverse of
558 // mean square distance (MSD)
559 for (lRefIndex=0; lRefIndex<lNumRefSamples; ++lRefIndex)
560 {
561 // Get comparison metric between the two samples
562 if ((enErrorNumber=ComparePhaseSamples(fpTargetSamples[lTargetIndex], fppRefSamples[lRefIndex],
563 irFrequencyRange, &fComparisonMetric))!=ERROR_NONE)
564 {
565 for (lRefIndex=0; lRefIndex<lNumRefSamples; ++lRefIndex) free(fppRefSamples[lRefIndex]);
566 free(fppRefSamples);
567 for (lTargetIndex=0; lTargetIndex<lNumRefSamples; ++lTargetIndex) free(fpTargetSamples[lTargetIndex]);
568 fclose(fpGraphFile);
569 free(fpTargetSamples);
570 free(graphFileName);
571 return enErrorNumber;
572 }
573
574 // Add to top 10 peaks if appropriate
575 sPeakIndex=10;
576 for (i=9; i>=0; --i)
577 {
578 if (fComparisonMetric<=faPeakValues[i]) break;
579 sPeakIndex=i;
580 }
581 if (sPeakIndex<10)
582 {
583 for (j=9; j>sPeakIndex; --j)
584 {
585 faPeakValues[j]=faPeakValues[j-1];
586 fvaPeakCoords[j].x=fvaPeakCoords[j-1].x;
587 fvaPeakCoords[j].y=fvaPeakCoords[j-1].y;
588 }
589 faPeakValues[sPeakIndex]=fComparisonMetric;
590 fvaPeakCoords[sPeakIndex].x=fIncrementInDegrees*lRefIndex;
591 fvaPeakCoords[sPeakIndex].y=fIncrementInDegrees*lTargetIndex;
592 }
593
594 // Write graph file entry
595 fprintf(fpGraphFile, ",%3.2f", fComparisonMetric);
596 }
597
598 // New line in graph file
599 fprintf(fpGraphFile, "\n");
600 }
601
602 // Close graph file
603 fclose(fpGraphFile);
604
605 // Free target samples
606 for (lTargetIndex=0; lTargetIndex<lNumRefSamples; ++lTargetIndex) free(fpTargetSamples[lTargetIndex]);
607 free(fpTargetSamples);
608
609 // Open analysis file
610 sprintf(csAnalysisFileName,"%s/%s_LOI_summ%ld.csv", csSearchPath, csRootName, lPlaneIndex);
611 if (bAnalysisFile && !(fpAnalysisFile=fopen(csAnalysisFileName, "w")))
612 return ErrorOpeningFile(csAnalysisFileName);
613
614 // Output analysis to analysis file
615 if (bAnalysisFile)
616 {
617 fprintf(fpAnalysisFile, "Rank\t1/MSE\tRef. Coord.\tTarget Coord.\n");
618 for (i=0, j=1; i<10; ++i)
619 {
620 fprintf(fpAnalysisFile, "%d\t%3.2f\t%3.2f\t%3.2f\n", j++,
621 faPeakValues[i], fvaPeakCoords[i].x, fvaPeakCoords[i].y);
622 }
623 }
624 *fpCoplanarAngularShift=fvaPeakCoords[0].y-fvaPeakCoords[0].x;
625 if (*fpCoplanarAngularShift>180) *fpCoplanarAngularShift-=360;
626 else if (*fpCoplanarAngularShift<-180) *fpCoplanarAngularShift+=360;
627
628 // Close analysis file
629 if (bAnalysisFile) fclose(fpAnalysisFile);
630 }
631
632 // Cleanup
633 for (lRefIndex=0; lRefIndex<lNumRefSamples; ++lRefIndex)
634 free(fppRefSamples[lRefIndex]);
635 free(fppRefSamples);
636 free(graphFileName);
637
638 /* TEMPORARAILY COMMENTED OUT 2020-09-28
639 // Error analysis
640 dCommonLinesAngleErr=fabs((double)(*fpCoplanarAngularShift-ExpectedAngularShift(Z_AXIS)));
641 while (dCommonLinesAngleErr>180) dCommonLinesAngleErr=360-dCommonLinesAngleErr;
642 while (dCommonLinesAngleErr<-180) dCommonLinesAngleErr=360+dCommonLinesAngleErr;
643 */
644
645 fCommonLinesPeak=faPeakValues[0];
646 for (i=0; i<10; ++i)
647 {
648 if (fabs((double)(*fpCoplanarAngularShift-GetAngularShift(fvaPeakCoords[i])))>2.0f)
649 {
650 fCommonLinesPeakDiff=faPeakValues[0]-faPeakValues[i];
651 break;
652 }
653 }
654
655 return ERROR_NONE;
656 }
657
658 /*********************************************************************/
CheckeredInversion(COMPLEX ** input,unsigned int dimension)659 COMPLEX **CheckeredInversion(COMPLEX **input, unsigned int dimension)
660 /*********************************************************************/
661 {
662 unsigned int i, j;
663
664 for (i=0; i<dimension; i+=2)
665 {
666 for (j=1; j<dimension; j+=2){
667 (*(input+i)+j)->imag = 0 - (*(input+i)+j)->imag;
668 (*(input+i)+j)->real = 0 - (*(input+i)+j)->real;
669 }
670 }
671
672 for (i=1; i<dimension; i+=2)
673 {
674 for (j=0; j<dimension; j+=2){
675 (*(input+i)+j)->imag = 0 - (*(input+i)+j)->imag;
676 (*(input+i)+j)->real = 0 - (*(input+i)+j)->real;
677 }
678 }
679
680 return(input);
681 }
682
outputRealAndImaginaryPlanes(COMPLEX ** TwoDFts,THD_3dim_dataset * din)683 int outputRealAndImaginaryPlanes(COMPLEX **TwoDFts, THD_3dim_dataset *din){
684 int x, y, inc;
685 char *prefix=DSET_PREFIX(din);
686 char *searchPath=DSET_DIRNAME(din);
687 char *outputFileName;
688 float *outData;
689 char appendageR[] = "Real", appendageI[] = "Imaginary";
690 COMPLEX *row;
691
692 // Determine input dimensions
693 int ny = DSET_NY(din);
694 int nx = DSET_NX(din);
695
696 // Make output array
697 if (!(outData=(float*)malloc(nx*ny*sizeof(float)))) return 0;
698
699 // Allocate memory to output name buffer
700 if (!(outputFileName=(char *)malloc(strlen(searchPath)+strlen(prefix)+64))){
701 free(outData);
702 return 0;
703 }
704
705 // Fill output data with real data
706 for (y=inc=0; y<ny; ++y){
707 row = TwoDFts[y];
708 for (x=0; x<nx; ++x)
709 outData[inc++]=row[x].real;
710 }
711
712 // Determine whether real file already exists
713 int outputFileExists = doesFileExist(searchPath,prefix,appendageR,outputFileName);
714
715 // Output Fourier spectrum image (if it does not already exist)
716 if (!outputFileExists){
717 THD_3dim_dataset *dout = EDIT_empty_copy(din);
718 sprintf(outputFileName,"%s%s%s",searchPath,prefix,appendageR);
719 EDIT_dset_items( dout ,
720 ADN_prefix, outputFileName,
721 ADN_none ) ;
722 EDIT_substitute_brick(dout, 0, MRI_float, outData);
723 DSET_write(dout);
724 }
725
726 // Fill output data with real data
727 for (y=inc=0; y<ny; ++y){
728 row = TwoDFts[y];
729 for (x=0; x<nx; ++x)
730 outData[inc++]=row[x].imag;
731 }
732
733 // Determine whether real file already exists
734 outputFileExists = doesFileExist(searchPath,prefix,appendageI,outputFileName);
735
736 // Output Fourier spectrum image (if it does not already exist)
737 if (!outputFileExists){
738 THD_3dim_dataset *dout = EDIT_empty_copy(din);
739 sprintf(outputFileName,"%s%s%s",searchPath,prefix,appendageI);
740 EDIT_dset_items( dout ,
741 ADN_prefix, outputFileName,
742 ADN_none ) ;
743 EDIT_substitute_brick(dout, 0, MRI_float, outData);
744 DSET_write(dout);
745 }
746
747 // Cleanup
748 free(outputFileName);
749 free(outData);
750
751 return 1;
752
753 }
754
outputFourierSpectrum(COMPLEX ** TwoDFts,THD_3dim_dataset * din)755 int outputFourierSpectrum(COMPLEX **TwoDFts, THD_3dim_dataset *din){
756 int x, y, inc;
757 char *prefix=DSET_PREFIX(din);
758 char *searchPath=DSET_DIRNAME(din);
759 char *outputFileName;
760 float *outData;
761 char appendage[] = "FTSpect";
762 COMPLEX *row;
763
764 // Determine input dimensions
765 int ny = DSET_NY(din);
766 int nx = DSET_NX(din);
767
768 // Make output array
769 if (!(outData=(float*)malloc(nx*ny*sizeof(float)))) return 0;
770
771 // Fill output data with Fourier spectral data
772 for (y=inc=0; y<ny; ++y){
773 row = TwoDFts[y];
774 for (x=0; x<nx; ++x)
775 outData[inc++]=hypot(row[x].real,row[x].imag);
776 }
777
778 // Allocate memory to output name buffer
779 if (!(outputFileName=(char *)malloc(strlen(searchPath)+strlen(prefix)+64))){
780 free(outData);
781 return 0;
782 }
783
784 // Determine whether output file already exists
785 int outputFileExists = doesFileExist(searchPath,prefix,appendage,outputFileName);
786
787 // Output Fourier spectrum image (if it does not already exist)
788 if (!outputFileExists){
789 THD_3dim_dataset *dout = EDIT_empty_copy(din);
790 sprintf(outputFileName,"%s%s%s",searchPath,prefix,appendage);
791 EDIT_dset_items( dout ,
792 ADN_prefix, outputFileName,
793 ADN_none ) ;
794 EDIT_substitute_brick(dout, 0, MRI_float, outData);
795 DSET_write(dout);
796 }
797
798 // Cleanup
799 free(outputFileName);
800 free(outData);
801
802 return 1;
803 }
804
get2DFourierTransform(THD_3dim_dataset * din,COMPLEX *** TwoDFts)805 int get2DFourierTransform(THD_3dim_dataset *din, COMPLEX ***TwoDFts){
806 int x, y, X, Y, inc;
807 float *inputImage;
808 COMPLEX *row, *newRow, **centeredFT;
809
810 // Get dimensions
811 int ny = DSET_NY(din);
812 int nx = DSET_NX(din);
813
814 // Allocate memory to complex Fourier plane
815 if (!(*TwoDFts=(COMPLEX **)malloc(ny*sizeof(COMPLEX *)))) return 0;
816 for (y=0; y<ny; ++y) if (!((*TwoDFts)[y]=(COMPLEX *)calloc(nx,sizeof(COMPLEX)))){
817 for (--y; y>=0; --y) free((*TwoDFts)[y]);
818 free((*TwoDFts));
819 return 0;
820 }
821
822 // Fill real components with spatial image data
823 inputImage = DSET_ARRAY(din, 0);
824 for (y=inc=0; y<ny; ++y){
825 row = (*TwoDFts)[y];
826 for (x=0; x<nx; ++x)
827 row[x].real = inputImage[inc++];
828 }
829
830 // Fouier transform plane
831 if (!FFT2D(*TwoDFts, nx, ny, 1)){
832 for (y=0; y<ny; ++y) free((*TwoDFts)[y]);
833 free(*TwoDFts);
834 *TwoDFts = NULL;
835 return 0;
836 }
837
838 // Allocate memory to centering buffer
839 if (!(centeredFT=(COMPLEX **)malloc(ny*sizeof(COMPLEX *)))) {
840 for (y=0; y<ny; ++y) free((*TwoDFts)[y]);
841 free(*TwoDFts);
842 *TwoDFts = NULL;
843 return 0;
844 }
845 for (y=0; y<ny; ++y) if (!(centeredFT[y]=(COMPLEX *)calloc(nx,sizeof(COMPLEX)))){
846 for (--y; y>=0; --y) free(centeredFT[y]);
847 free(centeredFT);
848 for (y=0; y<ny; ++y) free((*TwoDFts)[y]);
849 free(*TwoDFts);
850 *TwoDFts = NULL;
851 return 0;
852 }
853
854 // Center Fourier transform
855 for (y=0, Y=ny/2; Y<ny; ++y, ++Y){
856 // New UL and UR quadrant
857 newRow=centeredFT[Y];
858 row = (*TwoDFts)[y];
859 for (x=0, X=nx/2; X<nx; ++x, ++X){
860 newRow[x]=row[X]; // UL quadrant
861 newRow[X]=row[x]; // UR quadrant
862 }
863
864 // New LL and LR quadrant
865 newRow=centeredFT[y];
866 row = (*TwoDFts)[Y];
867 for (x=0, X=nx/2; X<nx; ++x, ++X){
868 newRow[x]=row[X]; // LL quadrant
869 newRow[X]=row[x]; // LR quadrant
870 }
871 }
872
873 // Free old, uncentered FT
874 for (y=0; y<ny; ++y) free((*TwoDFts)[y]);
875 free(*TwoDFts);
876
877 // Centered FT becomes new FT
878 *TwoDFts = centeredFT;
879
880 // Checkerboard invert real and imaginary planes
881 *TwoDFts=CheckeredInversion(*TwoDFts, ny);
882
883 return 1;
884 }
885
zeroPadProjection(THD_3dim_dataset * din,THD_3dim_dataset ** dout,int largestDimension,int factor)886 int zeroPadProjection(THD_3dim_dataset *din, THD_3dim_dataset **dout, int largestDimension, int factor){
887 char *prefix=DSET_PREFIX(din);
888 char *searchPath=DSET_DIRNAME(din);
889 char *outputFileName;
890 float *outData;
891 int outputDimension=largestDimension*factor;
892 int y, inInc;
893 char appendage[8];
894 THD_ivec3 iv_nxyz;
895
896 // Set output name appendage
897 sprintf(appendage, "pad%d", factor);
898
899 // Raise output dimension to the next integral power of two (for Fourier transform)
900 outputDimension = pow(2, ceil(log(outputDimension)/log(2)));
901
902 // Make output array
903 if (!(outData=(float*)calloc(outputDimension*outputDimension,sizeof(float)))) return 0;
904
905 // Determine input dimensions
906 int ny = DSET_NY(din);
907 int nx = DSET_NX(din);
908
909 // Determine where input begins in output
910 int xOffset=(int)((outputDimension-nx)/2);
911 int yOffset=(int)((outputDimension-ny)/2);
912 int offset=(outputDimension*yOffset)+xOffset;
913
914 // Get input array
915 float *indata = DSET_ARRAY(din, 0);
916
917 // Write input to output
918 size_t bytes2Copy = nx*sizeof(float);
919 for (y=inInc=0; y<ny; ++y){
920 memcpy((void *)&(outData[offset]), (void *)&(indata[inInc]), bytes2Copy);
921 offset+=outputDimension;
922 inInc+=nx;
923 }
924
925 // Allocate memory to output name buffer
926 if (!(outputFileName=(char *)malloc(strlen(searchPath)+strlen(prefix)+64))){
927 free(outData);
928 return 0;
929 }
930
931 // Determine whether output file already exists
932 int outputFileExists = doesFileExist(searchPath,prefix,appendage,outputFileName);
933
934 // Make output dataset
935 *dout = EDIT_empty_copy(din);
936 LOAD_IVEC3( iv_nxyz , outputDimension , outputDimension , 1 ) ;
937 sprintf(outputFileName,"%s%s%s",searchPath,prefix,appendage);
938 EDIT_dset_items( *dout ,
939 ADN_prefix, outputFileName,
940 ADN_nxyz , iv_nxyz ,
941 ADN_none ) ;
942 EDIT_substitute_brick(*dout, 0, MRI_float, outData);
943
944 // Output padded image to file
945 if( !outputFileExists ) { // If even file does not already exist
946 DSET_write(*dout);
947 }
948
949 // Cleanup
950 free(outputFileName);
951
952 return 1;
953 }
954
doesFileExist(char * searchPath,char * prefix,char * appendage,char * outputFileName)955 int doesFileExist(char * searchPath, char * prefix,char *appendage , char * outputFileName){
956 int outputFileExists=0;
957
958 struct dirent *dir;
959 DIR *d;
960 d = opendir(searchPath);
961 sprintf(outputFileName,"%s%s",prefix,appendage);
962 if (d) {
963 while ((dir = readdir(d)) != NULL) {
964 if (strstr(dir->d_name,outputFileName)){
965 outputFileExists=1;
966 break;
967 }
968 }
969 closedir(d);
970 }
971
972 return outputFileExists;
973 }
974
float2DImage(THD_3dim_dataset * dset)975 int float2DImage(THD_3dim_dataset *dset){
976 char *prefix=DSET_PREFIX(dset);
977 char *searchPath=DSET_DIRNAME(dset);
978 char *outputFileName;
979 char appendage[8];
980
981 // Subtract mean of perimeter from entire image to avoid spurious edge effects
982 // in Fourier transform
983 float mean=0;
984 int count=0, x, y;
985
986 // Get input pixel data
987 float *indata = DSET_ARRAY(dset, 0);
988
989 // Determine input dimensions
990 int ny = DSET_NY(dset);
991 int nx = DSET_NX(dset);
992 int lastY = ny-1;
993 int lastX = nx-1;
994
995 // Get mean value arount perimeter
996 int offset=lastY*nx;
997 for (x=0; x<nx; ++x){
998 mean += indata[x]+indata[offset++];
999 count+=2;
1000 }
1001 offset=nx;
1002 for (y=0; y<lastY; ++y){
1003 mean += indata[offset];
1004 offset += lastX;
1005 mean += indata[offset++];
1006 count+=2;
1007 }
1008 mean /= count;
1009
1010 // Subtract perimeter mean from whole image
1011 int pixelCount=nx*ny;
1012 for (x=0; x<pixelCount; ++x) indata[x]-=mean;
1013
1014 // Allocate memory to output name buffer
1015 if (!(outputFileName=(char *)malloc(strlen(searchPath)+strlen(prefix)+64))){
1016 return 0;
1017 }
1018
1019 // Set output name appendage
1020 sprintf(appendage, "float");
1021
1022 // See if image exists
1023 int outputFileExists = doesFileExist(searchPath,prefix,appendage,outputFileName);
1024
1025 // Output floated image to file
1026 sprintf(outputFileName,"%s%s%s",searchPath,prefix,appendage);
1027 EDIT_dset_items( dset ,
1028 ADN_prefix, outputFileName,
1029 ADN_none ) ;
1030 if( !outputFileExists ) { // If even file does not already exist
1031 DSET_write(dset);
1032 }
1033
1034 // Cleanup
1035 free(outputFileName);
1036
1037 return 1;
1038 }
1039
getLargestDimension(THD_3dim_dataset * din)1040 int getLargestDimension(THD_3dim_dataset *din){
1041
1042 // Determine input dimensions
1043 int nz = DSET_NZ(din);
1044 int ny = DSET_NY(din);
1045 int nx = DSET_NX(din);
1046
1047 return (nx>ny)? ((nx>nz)? nx : nz) : ((ny>nz)? ny : nz);
1048 }
1049
makeProjection(THD_3dim_dataset * din,THD_3dim_dataset ** dout,char projCode,char * orientation)1050 int makeProjection(THD_3dim_dataset *din, THD_3dim_dataset **dout, char projCode, char *orientation){
1051 int nz, ny, nx, inInc, outInc, rows, cols, x;
1052 char *prefix=DSET_PREFIX(din);
1053 char *searchPath=DSET_DIRNAME(din);
1054 char *outputCode=(projCode=='a')? "Axial" : ((projCode=='c')? "Coronal" : "Sagital");
1055 char *outputFileName;
1056 int outputFileExists=0, outPixelCount, start;
1057 THD_ivec3 iv_nxyz;
1058 float *outData=NULL;
1059
1060 // Determine input dimensions
1061 nz = DSET_NZ(din);
1062 ny = DSET_NY(din);
1063 nx = DSET_NX(din);
1064
1065 // Get input voxel data
1066 float *indata = DSET_ARRAY(din, 0);
1067
1068 // Determine axis
1069 switch (tolower(projCode)){
1070 case 'a':
1071 if (orientation[0]=='S' || orientation[0]=='I') projCode='x';
1072 else if (orientation[1]=='S' || orientation[1]=='I') projCode='y';
1073 else if (orientation[2]=='S' || orientation[2]=='I') projCode='z';
1074 else return ERROR_SEARCH;
1075 break;
1076 case 'c':
1077 if (orientation[0]=='A' || orientation[0]=='P') projCode='x';
1078 else if (orientation[1]=='A' || orientation[1]=='P') projCode='y';
1079 else if (orientation[2]=='A' || orientation[2]=='P') projCode='z';
1080 else return ERROR_SEARCH;
1081 break;
1082 case 's':
1083 if (orientation[0]=='L' || orientation[0]=='R') projCode='x';
1084 else if (orientation[1]=='L' || orientation[1]=='R') projCode='y';
1085 else if (orientation[2]=='L' || orientation[2]=='R') projCode='z';
1086 else return ERROR_SEARCH;
1087 break;
1088 }
1089
1090 // Apply required projection
1091 switch (projCode){
1092 case 'z': // Axial projection
1093
1094 // Determine output dimansions
1095 rows=ny;
1096 cols=nx;
1097
1098 // Allocate memory to output voxel data
1099 outPixelCount=rows*cols;
1100 if (!(outData=(float *)calloc(outPixelCount,sizeof(float)))) return 0;
1101
1102 // Fill output data
1103 start=0;
1104 for (inInc=0, outInc=0; outInc<outPixelCount; ++outInc){
1105 inInc=start;
1106 for (x=0; x<nz; ++x){
1107 outData[outInc]+=indata[inInc];
1108 inInc+=outPixelCount;
1109 }
1110 start+=1;
1111 }
1112 break;
1113 case 'y': // Coronal projection
1114
1115 // Determine output dimansions
1116 rows=nz;
1117 cols=nx;
1118
1119 // Allocate memory to output voxel data
1120 outPixelCount=rows*cols;
1121 if (!(outData=(float *)calloc(outPixelCount,sizeof(float)))) return 0;
1122
1123 // Fill output data
1124 start=0;
1125 int z, i;
1126 for (outInc=outPixelCount-1, z=0; z<nz; ++z){
1127 start=z*nx*ny;
1128 for (x=0; x<nx; ++x){
1129 inInc=start;
1130 for (i=0; i<ny; ++i){
1131 outData[outInc]+=indata[inInc];
1132 inInc+=nx;
1133 }
1134 start+=1;
1135 --outInc;
1136 }
1137 }
1138 break;
1139 case 'x': // Sagital projection
1140
1141 // Determine output dimansions
1142 rows=nz;
1143 cols=ny;
1144
1145 // Allocate memory to output voxel data
1146 outPixelCount=rows*cols;
1147 if (!(outData=(float *)calloc(outPixelCount,sizeof(float)))) return 0;
1148
1149 // Fill output data
1150 for (inInc=0, outInc=rows*cols-1; outInc>=0; --outInc){
1151 for (x=0; x<nx; ++x) outData[outInc]+=indata[inInc++];
1152 }
1153 break;
1154 }
1155
1156 // Allocate memory to output name buffer
1157 if (!(outputFileName=(char *)malloc(strlen(searchPath)+strlen(prefix)+64))){
1158 free(outData);
1159 return 0;
1160 }
1161
1162 // Determine whether output file already exists
1163 struct dirent *dir;
1164 DIR *d;
1165 d = opendir(searchPath);
1166 sprintf(outputFileName,"%s%s",prefix,outputCode);
1167 if (d) {
1168 while ((dir = readdir(d)) != NULL) {
1169 if (strstr(dir->d_name,outputFileName)){
1170 outputFileExists=1;
1171 }
1172 }
1173 closedir(d);
1174 }
1175
1176 /* Output 2D projection volume */
1177 *dout = EDIT_empty_copy(din);
1178 LOAD_IVEC3( iv_nxyz , cols , rows , 1 ) ;
1179 sprintf(outputFileName,"%s%s%s",searchPath,prefix,outputCode);
1180 EDIT_dset_items( *dout ,
1181 ADN_prefix, outputFileName,
1182 ADN_nxyz , iv_nxyz ,
1183 ADN_none ) ;
1184 EDIT_substitute_brick(*dout, 0, MRI_float, outData);
1185 if( !outputFileExists ) { // If even file does not already exist
1186 DSET_write(*dout);
1187 }
1188
1189 // Cleanup (Don't free outData)
1190 free(outputFileName);
1191
1192 return 1;
1193 }
1194
shortToFloat(THD_3dim_dataset ** din)1195 ERROR_NUMBER shortToFloat(THD_3dim_dataset **din){
1196 float *outputBuffer;
1197 short *indata;
1198 int i;
1199
1200 // Determine input dimensions
1201 int nz = DSET_NZ(*din);
1202 int ny = DSET_NY(*din);
1203 int nx = DSET_NX(*din);
1204 int numVoxels=nx*ny*nz;
1205
1206 // Memory allocation
1207 if (!(outputBuffer=(float *)malloc(numVoxels*sizeof(float)))){
1208 return ERROR_MEMORY_ALLOCATION;
1209 }
1210
1211 // Fill input data
1212 indata = DSET_ARRAY(*din, 0);
1213
1214 // Convert short int data to floating point
1215 for (i=0; i<numVoxels; ++i){
1216 outputBuffer[i]=indata[i];
1217 }
1218
1219 // Update dataset
1220 EDIT_dset_items( *din ,
1221 ADN_type, MRI_float,
1222 ADN_none ) ;
1223 EDIT_substitute_brick(*din, 0, MRI_float, outputBuffer);
1224
1225 return ERROR_NONE;
1226 }
1227
unweave_sections(THD_3dim_dataset * din,THD_3dim_dataset ** dodd,THD_3dim_dataset ** deven)1228 int unweave_sections(THD_3dim_dataset *din, THD_3dim_dataset **dodd, THD_3dim_dataset **deven){
1229 char *outputFileName;
1230 int i, nz, ny, nx, nodd, neven, readInc, inIndex, outIndex, planeSize;
1231 size_t bytes2Copy;
1232 float * indata, *oddData=NULL, *evenData=NULL;
1233 THD_ivec3 iv_nxyz;
1234 int bEvenFileExists=0, bOddFileExists=0;
1235
1236 // Determine input dimensions
1237 nz = DSET_NZ(din);
1238 ny = DSET_NY(din);
1239 nx = DSET_NX(din);
1240
1241 // Determine output section counts
1242 neven=(int)((nz+1)/2);
1243 nodd=nz-neven;
1244
1245 // Determine plane size read increment
1246 planeSize = nx*ny;
1247 readInc = 2*planeSize;
1248
1249 // Allocate memory to output voxel data
1250 bytes2Copy=planeSize*sizeof(float);
1251 if (!(oddData=(float *)malloc(nodd*bytes2Copy)) ||
1252 !(evenData=(float *)malloc(neven*bytes2Copy))){
1253 if (oddData) free(oddData);
1254 return 0;
1255 }
1256
1257 // Fill even data
1258 indata = DSET_ARRAY(din, 0);
1259 inIndex=outIndex=0;
1260 for (i=0; i<neven; ++i){
1261 memcpy((void *)&evenData[outIndex],(void *)&indata[inIndex],bytes2Copy);
1262 inIndex+=readInc;
1263 outIndex+=planeSize;
1264 }
1265
1266 // Fill odd data
1267 outIndex=0;
1268 inIndex=planeSize;
1269 for (i=0; i<nodd; ++i){
1270 memcpy((void *)&oddData[outIndex],(void *)&indata[inIndex],bytes2Copy);
1271 inIndex+=readInc;
1272 outIndex+=planeSize;
1273 }
1274
1275 // Get search path, and prefix, from dataset
1276 char *prefix=DSET_PREFIX(din);
1277 char *searchPath=DSET_DIRNAME(din);
1278
1279 // Allocate memory to output name buffer
1280 if (!(outputFileName=(char *)malloc(strlen(searchPath)+strlen(prefix)+32))){
1281 free(oddData);
1282 free(evenData);
1283 return 0;
1284 }
1285
1286 // Determine whether files actually exist
1287 struct dirent *dir;
1288 DIR *d;
1289 d = opendir(searchPath);
1290 if (d) {
1291 while ((dir = readdir(d)) != NULL) {
1292 sprintf(outputFileName,"%sEven",prefix);
1293 if (strstr(dir->d_name,outputFileName)){
1294 bEvenFileExists=1;
1295 } else {
1296 sprintf(outputFileName,"%sOdd",prefix);
1297 if (strstr(dir->d_name,outputFileName)){
1298 bOddFileExists=1;
1299 }
1300 }
1301 }
1302 closedir(d);
1303 }
1304
1305
1306 /* Output even volume */
1307 *deven = EDIT_empty_copy(din);
1308 LOAD_IVEC3( iv_nxyz , nx , ny , neven ) ;
1309 sprintf(outputFileName,"%s%sEven",searchPath,prefix);
1310 EDIT_dset_items( *deven ,
1311 ADN_prefix, outputFileName,
1312 ADN_nxyz , iv_nxyz ,
1313 ADN_none ) ;
1314 EDIT_substitute_brick(*deven, 0, MRI_float, evenData);
1315 if( !bEvenFileExists ) { // If even file does not already exist
1316 DSET_write(*deven);
1317 }
1318
1319 *dodd = EDIT_empty_copy(din);
1320 LOAD_IVEC3( iv_nxyz , nx , ny , nodd ) ;
1321 sprintf(outputFileName,"%s%sOdd",searchPath,prefix);
1322 EDIT_dset_items( *dodd ,
1323 ADN_prefix, outputFileName,
1324 ADN_nxyz , iv_nxyz ,
1325 ADN_none ) ;
1326 EDIT_substitute_brick(*dodd, 0, MRI_float, oddData);
1327 if( !bOddFileExists ) { // If odd file does not already exist
1328 DSET_write(*dodd);
1329 }
1330
1331 // Cleanup (Don't free evenData or oddData)
1332 free(outputFileName);
1333
1334 return 1;
1335 }
1336
open_input_dset(THD_3dim_dataset ** din,char * fname)1337 int open_input_dset(THD_3dim_dataset ** din, char * fname)
1338 {
1339 int errorNumber;
1340 *din = THD_open_dataset(fname);
1341 if( ! *din ) {
1342 fprintf(stderr,"** failed to read input dataset '%s'\n", fname);
1343 return ERROR_READING_FILE;
1344 }
1345
1346 /* refuse to work with anything but float here */
1347 // int brickType=DSET_BRICK_TYPE(*din, 0);
1348 int brickType=DSET_BRICK_TYPE(*din, 0);
1349 if( brickType == MRI_float ) {
1350 /* data is not automatically read in, do it now */
1351 DSET_load(*din);
1352 /*
1353 fprintf(stderr,"** input must be of type float, failing...\n");
1354 return 1;
1355 */
1356 } else if (brickType==MRI_short){
1357
1358 DSET_load(*din);
1359
1360 if ((errorNumber=shortToFloat(din))!=ERROR_NONE)
1361 return errorNumber;
1362 // TODO: Add code
1363 }
1364
1365 return ERROR_NONE;
1366 }
1367
Cleanup(char * inputFileName,COMPLEX *** TwoDFts,THD_3dim_dataset * din)1368 int Cleanup(char *inputFileName, COMPLEX ***TwoDFts, THD_3dim_dataset *din){
1369 int i, j, ny;
1370
1371 if (inputFileName) free(inputFileName);
1372
1373 for (i=0; i<6; ++i){
1374 if (TwoDFts[i]){
1375 ny = DSET_NY(din);
1376 for (j=0; j<ny; ++j) free(TwoDFts[i][j]);
1377 free(TwoDFts[i]);
1378 }
1379 }
1380
1381 return 0;
1382 }
1383
1384
1385 /************************************ 2D FFT by Paul Bourke, 1993 ***************************/
1386
1387
1388
1389 /*-------------------------------------------------------------------------
1390 Perform a 2D FFT inplace given a complex 2D array
1391 The direction dir, 1 for forward, -1 for reverse
1392 The size of the array (nx,ny)
1393 Return false if there are memory problems or
1394 the dimensions are not powers of 2
1395 */
FFT2D(COMPLEX ** c,int nx,int ny,int dir)1396 int FFT2D(COMPLEX **c,int nx,int ny,int dir)
1397 {
1398 int i,j;
1399 int m,twopm;
1400 double *real,*imag;
1401
1402 /* Transform the rows */
1403 real = (double *)malloc(nx * sizeof(double));
1404 imag = (double *)malloc(nx * sizeof(double));
1405 if (real == NULL || imag == NULL)
1406 return(FALSE);
1407 if (!Powerof2(nx,&m,&twopm) || twopm != nx)
1408 return(FALSE);
1409 for (j=0;j<ny;j++) {
1410 for (i=0;i<nx;i++) {
1411 real[i] = c[i][j].real;
1412 imag[i] = c[i][j].imag;
1413 }
1414 FFT(dir,m,real,imag);
1415 for (i=0;i<nx;i++) {
1416 c[i][j].real = real[i];
1417 c[i][j].imag = imag[i];
1418 }
1419 }
1420 free(real);
1421 free(imag);
1422
1423 /* Transform the columns */
1424 real = (double *)malloc(ny * sizeof(double));
1425 imag = (double *)malloc(ny * sizeof(double));
1426 if (real == NULL || imag == NULL)
1427 return(FALSE);
1428 if (!Powerof2(ny,&m,&twopm) || twopm != ny)
1429 return(FALSE);
1430 for (i=0;i<nx;i++) {
1431 for (j=0;j<ny;j++) {
1432 real[j] = c[i][j].real;
1433 imag[j] = c[i][j].imag;
1434 }
1435 FFT(dir,m,real,imag);
1436 for (j=0;j<ny;j++) {
1437 c[i][j].real = real[j];
1438 c[i][j].imag = imag[j];
1439 }
1440 }
1441 free(real);
1442 free(imag);
1443
1444 return(TRUE);
1445 }
1446
1447 /*-------------------------------------------------------------------------
1448 This computes an in-place complex-to-complex FFT
1449 x and y are the real and imaginary arrays of 2^m points.
1450 dir = 1 gives forward transform
1451 dir = -1 gives reverse transform
1452
1453 Formula: forward
1454 N-1
1455 ---
1456 1 \ - j k 2 pi n / N
1457 X(n) = --- > x(k) e = forward transform
1458 N / n=0..N-1
1459 ---
1460 k=0
1461
1462 Formula: reverse
1463 N-1
1464 ---
1465 \ j k 2 pi n / N
1466 X(n) = > x(k) e = forward transform
1467 / n=0..N-1
1468 ---
1469 k=0
1470 */
FFT(int dir,int m,double * x,double * y)1471 int FFT(int dir,int m,double *x,double *y)
1472 {
1473 long nn,i,i1,j,k,i2,l,l1,l2;
1474 double c1,c2,tx,ty,t1,t2,u1,u2,z;
1475
1476 /* Calculate the number of points */
1477 nn = 1;
1478 for (i=0;i<m;i++)
1479 nn *= 2;
1480
1481 /* Do the bit reversal */
1482 i2 = nn >> 1;
1483 j = 0;
1484 for (i=0;i<nn-1;i++) {
1485 if (i < j) {
1486 tx = x[i];
1487 ty = y[i];
1488 x[i] = x[j];
1489 y[i] = y[j];
1490 x[j] = tx;
1491 y[j] = ty;
1492 }
1493 k = i2;
1494 while (k <= j) {
1495 j -= k;
1496 k >>= 1;
1497 }
1498 j += k;
1499 }
1500
1501 /* Compute the FFT */
1502 c1 = -1.0;
1503 c2 = 0.0;
1504 l2 = 1;
1505 for (l=0;l<m;l++) {
1506 l1 = l2;
1507 l2 <<= 1;
1508 u1 = 1.0;
1509 u2 = 0.0;
1510 for (j=0;j<l1;j++) {
1511 for (i=j;i<nn;i+=l2) {
1512 i1 = i + l1;
1513 t1 = u1 * x[i1] - u2 * y[i1];
1514 t2 = u1 * y[i1] + u2 * x[i1];
1515 x[i1] = x[i] - t1;
1516 y[i1] = y[i] - t2;
1517 x[i] += t1;
1518 y[i] += t2;
1519 }
1520 z = u1 * c1 - u2 * c2;
1521 u2 = u1 * c2 + u2 * c1;
1522 u1 = z;
1523 }
1524 c2 = sqrt((1.0 - c1) / 2.0);
1525 if (dir == 1)
1526 c2 = -c2;
1527 c1 = sqrt((1.0 + c1) / 2.0);
1528 }
1529
1530 /* Scaling for forward transform */
1531 if (dir == 1) {
1532 for (i=0;i<nn;i++) {
1533 x[i] /= (double)nn;
1534 y[i] /= (double)nn;
1535 }
1536 }
1537
1538 return(TRUE);
1539 }
1540
1541 /*-------------------------------------------------------------------------
1542 Calculate the closest but lower power of two of a number
1543 twopm = 2**m <= n
1544 Return TRUE if 2**m == n
1545 */
Powerof2(int n,int * m,int * twopm)1546 int Powerof2(int n,int *m,int *twopm)
1547 {
1548 if (n <= 1) {
1549 *m = 0;
1550 *twopm = 1;
1551 return(FALSE);
1552 }
1553
1554 *m = 1;
1555 *twopm = 2;
1556 do {
1557 (*m)++;
1558 (*twopm) *= 2;
1559 } while (2*(*twopm) <= n);
1560
1561 if (*twopm != n)
1562 return(FALSE);
1563 else
1564 return(TRUE);
1565 }
1566
1567 /**********************************************************************/
MakeComplexPlane(long iRows,long iColumns,COMPLEX ** Data)1568 ComplexPlane MakeComplexPlane(long iRows, long iColumns, COMPLEX **Data)
1569 /**********************************************************************/
1570 {
1571 ComplexPlane buffer;
1572 size_t stBytesPerRow=iColumns * sizeof(COMPLEX);
1573
1574 buffer.iRows=iRows;
1575 buffer.iColumns=iColumns;
1576 buffer.Data=Data;
1577
1578 if (Data == NULL)
1579 {
1580 if ((buffer.Data = (COMPLEX **)malloc(buffer.iRows * sizeof(COMPLEX *))))
1581 {
1582 for (long i=0; i<buffer.iRows; ++i)
1583 if ((buffer.Data[i]=(COMPLEX *)calloc(1, stBytesPerRow)) == NULL)
1584 {
1585 for (--i; i>=0; --i) free(buffer.Data[i]);
1586 buffer.Data = NULL;
1587 return(buffer);
1588 }
1589 }
1590 }
1591
1592 return(buffer);
1593 }
1594
1595 /*************************************************/
FreeComplexPlane(ComplexPlane Complex_Plane)1596 void FreeComplexPlane(ComplexPlane Complex_Plane)
1597 /*************************************************/
1598 {
1599 int i;
1600
1601 for (i=0; i<Complex_Plane.iRows; ++i) free(Complex_Plane.Data[i]);
1602 free(Complex_Plane.Data);
1603 Complex_Plane.Data = NULL;
1604 }
1605
1606 /**************************************************************************/
CopyComplexPlane(ComplexPlane cvInput,ComplexPlane cvOutput)1607 ERROR_NUMBER CopyComplexPlane(ComplexPlane cvInput, ComplexPlane cvOutput)
1608 /**************************************************************************/
1609 {
1610 int iBytes2Copy, i;
1611
1612 if ((cvOutput.iRows != cvInput.iRows) || (cvOutput.iColumns != cvInput.iColumns)) return ERROR_IMAGE_DIMENSIONS;
1613
1614 iBytes2Copy = cvInput.iColumns * sizeof(COMPLEX);
1615 for (i=0; i<cvInput.iRows; ++i)
1616 memcpy( (void *)(cvOutput.Data[i]), (void *)(cvInput.Data[i]), iBytes2Copy );
1617
1618 return ERROR_NONE;
1619 }
1620
RotScaleComplexPlane(ComplexPlane * cppComplexPlane,FloatVector fvTranslation,float fClockwiseRotationInDegrees,float fScaleX,float fScaleY,BOOL bClip,BYTE bInterpOrder)1621 ERROR_NUMBER RotScaleComplexPlane(ComplexPlane *cppComplexPlane, FloatVector fvTranslation,
1622 float fClockwiseRotationInDegrees, float fScaleX, float fScaleY, BOOL bClip, BYTE bInterpOrder)
1623 {
1624 Matrix mTransMatrix;
1625 float faaMatrixX[9];
1626 long lMatrixIndex=0;
1627 double dRotationInRadians=DEGREES2RADIANS*fClockwiseRotationInDegrees;
1628 // Reverse rotation since matrix is supposed to give source of data, not destination.
1629 double dCosTheta=cos(-dRotationInRadians), dSinTheta=sin(-dRotationInRadians);
1630
1631 // If image not being clipped, determine output dimensions
1632 if (!bClip) DetermineUnclippedOutputDimensions(cppComplexPlane->iColumns, cppComplexPlane->iRows,
1633 fClockwiseRotationInDegrees, fvTranslation, 1.0f/fScaleX, 1.0f/fScaleY);
1634
1635 // Reverse translation since matrix is supposed to give source of data, not destination.
1636 fvTranslation.x = -fvTranslation.x;
1637 fvTranslation.y = -fvTranslation.y;
1638
1639 faaMatrixX[lMatrixIndex++] = fScaleX*(float)dCosTheta;
1640 faaMatrixX[lMatrixIndex++] = -fScaleX*(float)dSinTheta;
1641 faaMatrixX[lMatrixIndex++] = fScaleX*fvTranslation.x;
1642 faaMatrixX[lMatrixIndex++] = fScaleY*(float)dSinTheta;
1643 faaMatrixX[lMatrixIndex++] = fScaleY*(float)dCosTheta;
1644 faaMatrixX[lMatrixIndex++] = fScaleY*fvTranslation.y;
1645 faaMatrixX[lMatrixIndex++] = 0;
1646 faaMatrixX[lMatrixIndex++] = 0;
1647 faaMatrixX[lMatrixIndex++] = 1;
1648 mTransMatrix.fpData=&(faaMatrixX[0]);
1649 mTransMatrix.lRows=mTransMatrix.lColumns=3;
1650
1651 return RotScaleComplexPlaneM(cppComplexPlane, mTransMatrix, bClip, bInterpOrder);
1652 }
1653
RotScaleComplexPlaneM(ComplexPlane * cppComplexPlane,Matrix mTransMatrix,BOOL bClip,BYTE bInterpOrder)1654 ERROR_NUMBER RotScaleComplexPlaneM(ComplexPlane *cppComplexPlane, Matrix mTransMatrix, BOOL bClip, BYTE bInterpOrder)
1655 {
1656 long i, j;
1657 long iMinX=bInterpOrder+1, iMaxX=cppComplexPlane->iColumns-bInterpOrder-1;
1658 long iMinY=bInterpOrder+1, iMaxY=cppComplexPlane->iRows-bInterpOrder-1;
1659 float **fppAmplitudeBuffer, **fppPhaseBuffer, *fpAmplitudeRow, *fpPhaseRow;
1660 float **fppOutputAmplitudeBuffer, **fppOutputPhaseBuffer, *fpOutputAmplitudeRow, *fpOutputPhaseRow;
1661 float fpInputVector[3], fpOutputVector[3];
1662 long lHalfRows=cppComplexPlane->iRows/2, lHalfCols=cppComplexPlane->iColumns/2;
1663 double dReal, dImaginary;
1664 COMPLEX *cpOutputRow;
1665 LongVector lvOriginalDimensions={cppComplexPlane->iColumns, cppComplexPlane->iRows};
1666
1667 // Memory must have been allocated to input plane
1668 if (!(cppComplexPlane->Data)) return ERROR_NULL_PTR;
1669
1670 // If clipping, output dimensions same as those of input
1671 if (bClip)
1672 {
1673 lvOutputDimensions.x=cppComplexPlane->iRows;
1674 lvOutputDimensions.y=cppComplexPlane->iColumns;
1675 }
1676
1677 // Allocate memory to phase and amplitude buffers
1678 if (!(fppAmplitudeBuffer=(float **)malloc(cppComplexPlane->iRows*sizeof(float *))) ||
1679 !(fppPhaseBuffer=(float **)malloc(cppComplexPlane->iRows*sizeof(float *))) ||
1680 !(fppOutputAmplitudeBuffer=(float **)malloc(lvOutputDimensions.y*sizeof(float *))) ||
1681 !(fppOutputPhaseBuffer=(float **)malloc(lvOutputDimensions.y*sizeof(float *))))
1682 {
1683 if (fppAmplitudeBuffer) free(fppAmplitudeBuffer);
1684 if (fppPhaseBuffer) free(fppPhaseBuffer);
1685 if (fppOutputAmplitudeBuffer) free(fppOutputAmplitudeBuffer);
1686 return ERROR_MEMORY_ALLOCATION;
1687 }
1688 for (i=0; i<cppComplexPlane->iRows; ++i)
1689 if (!(fppAmplitudeBuffer[i]=(float *)malloc(cppComplexPlane->iColumns*sizeof(float))) ||
1690 !(fppPhaseBuffer[i]=(float *)malloc(cppComplexPlane->iColumns*sizeof(float))))
1691 {
1692 if (fppAmplitudeBuffer[i]) free(fppAmplitudeBuffer[i]);
1693 for (--i; i>=0; --i)
1694 {
1695 free(fppAmplitudeBuffer[i]);
1696 free(fppPhaseBuffer[i]);
1697 }
1698 free(fppAmplitudeBuffer);
1699 free(fppPhaseBuffer);
1700 free(fppOutputAmplitudeBuffer);
1701 free(fppOutputPhaseBuffer);
1702 return ERROR_MEMORY_ALLOCATION;
1703 }
1704 for (i=0; i<lvOutputDimensions.y; ++i)
1705 if (!(fppOutputAmplitudeBuffer[i]=(float *)malloc(lvOutputDimensions.x*sizeof(float))) ||
1706 !(fppOutputPhaseBuffer[i]=(float *)malloc(lvOutputDimensions.x*sizeof(float))))
1707 {
1708 if (fppOutputAmplitudeBuffer[i]) free(fppOutputAmplitudeBuffer[i]);
1709 for (--i; i>=0; --i)
1710 {
1711 free(fppOutputAmplitudeBuffer[i]);
1712 free(fppOutputPhaseBuffer[i]);
1713 }
1714 for (i=0; i<cppComplexPlane->iRows; ++i)
1715 {
1716 free(fppAmplitudeBuffer[i]);
1717 free(fppPhaseBuffer[i]);
1718 }
1719 free(fppAmplitudeBuffer);
1720 free(fppPhaseBuffer);
1721 free(fppOutputAmplitudeBuffer);
1722 free(fppOutputPhaseBuffer);
1723 return ERROR_MEMORY_ALLOCATION;
1724 }
1725
1726 // Load phase and amplitude buffers
1727 for (i=0; i<cppComplexPlane->iRows; ++i)
1728 {
1729 fpAmplitudeRow=fppAmplitudeBuffer[i];
1730 fpPhaseRow=fppPhaseBuffer[i];
1731 for (j=0; j<cppComplexPlane->iColumns; ++j)
1732 {
1733 dReal=cppComplexPlane->Data[i][j].real;
1734 dImaginary=cppComplexPlane->Data[i][j].imag;
1735 fpAmplitudeRow[j] = (float)Pythag(dReal, dImaginary);
1736 fpPhaseRow[j] = (float)((dImaginary==0)? 0 : ((dReal==0.0)? ((dImaginary>0)? MATH_PI_2 : -MATH_PI_2 ) : atan((dImaginary/dReal))));
1737 if (dReal<0)
1738 {
1739 fpPhaseRow[j]+=(float)MATH_PI;
1740 if (fpPhaseRow[j]>(float)CIRCLE) fpPhaseRow[j]-=(float)CIRCLE;
1741 else if (fpPhaseRow[j]<0) fpPhaseRow[j]+=(float)CIRCLE;
1742 }
1743 }
1744 }
1745
1746 if (bClip)
1747 {
1748 fpInputVector[2] = 1;
1749 for (i=0; i<cppComplexPlane->iRows; ++i)
1750 {
1751 fpOutputAmplitudeRow=fppOutputAmplitudeBuffer[i];
1752 fpOutputPhaseRow=fppOutputPhaseBuffer[i];
1753
1754 fpInputVector[1] = (float)(i-lHalfRows); // Rotation around centre of image
1755 for (j=0; j<cppComplexPlane->iColumns; ++j)
1756 {
1757 fpInputVector[0] = (float)(j-lHalfCols); // Rotation around centre of image
1758
1759 // Determine transformation from point on output image to corresponding point on input image
1760 MatVectMult(mTransMatrix, &(fpInputVector[0]), &(fpOutputVector[0]));
1761
1762 // Correct for offset from UL corner to centre
1763 fpOutputVector[0] += lHalfCols;
1764 fpOutputVector[1] += lHalfRows;
1765
1766 // Find corresponding value if input image point within input image
1767 if (bInterpOrder>NEAREST_NEIGHBOR_INTERPOLATION && fpOutputVector[0]>iMinX && fpOutputVector[0]<iMaxX &&
1768 fpOutputVector[1]>iMinY && fpOutputVector[1]<iMaxY)
1769 {
1770 fpOutputAmplitudeRow[j] = (
1771 (bInterpOrder==BILINEAR_INTERPOLATION)? (float)Interp_BL(fppAmplitudeBuffer, fpOutputVector[0], fpOutputVector[1]) :
1772 ((bInterpOrder==BIQUADRATIC_INTERPOLATION)? (float)Interp_BQ(fppAmplitudeBuffer, fpOutputVector[0], fpOutputVector[1]) :
1773 (float)Interp_BC(fppAmplitudeBuffer, fpOutputVector[0], fpOutputVector[1])));
1774 fpOutputPhaseRow[j] = (
1775 (bInterpOrder==BILINEAR_INTERPOLATION)? (float)Interp_BL(fppPhaseBuffer, fpOutputVector[0], fpOutputVector[1]) :
1776 ((bInterpOrder==BIQUADRATIC_INTERPOLATION)? (float)Interp_BQ(fppPhaseBuffer, fpOutputVector[0], fpOutputVector[1]) :
1777 (float)Interp_BC(fppPhaseBuffer, fpOutputVector[0], fpOutputVector[1])));
1778 }
1779 else
1780 {
1781 FixedPtVector iNearestNeighbor={Round(fpOutputVector[0]), Round(fpOutputVector[1])};
1782
1783 if (iNearestNeighbor.x>=0 && iNearestNeighbor.x<cppComplexPlane->iColumns &&
1784 iNearestNeighbor.y>=0 && iNearestNeighbor.y<cppComplexPlane->iRows)
1785 {
1786 fpOutputAmplitudeRow[j] = (fppAmplitudeBuffer[iNearestNeighbor.y][iNearestNeighbor.x]);
1787 fpOutputPhaseRow[j] = (fppPhaseBuffer[iNearestNeighbor.y][iNearestNeighbor.x]);
1788 }
1789 else
1790 fpOutputAmplitudeRow[j] = fpOutputPhaseRow[j] = 0;
1791 }
1792 }
1793 }
1794 }
1795 else
1796 {
1797 LongRange lrXRange, lrYRange;
1798 long lLargeX, lLargeY, lXOffset, lYOffset, x, y;
1799
1800 // Ensure output dimensions have been determined
1801 if (lvOutputDimensions.x==-1) return ERROR_IMAGE_DIMENSIONS;
1802
1803 // Determine dimensions of large rectangle that can just hold both the input and output images
1804 lLargeX=fmax(lvOutputDimensions.x, cppComplexPlane->iColumns);
1805 lLargeY=fmax(lvOutputDimensions.y, cppComplexPlane->iRows);
1806 lHalfCols=lLargeX/2;
1807 lHalfRows=lLargeY/2;
1808
1809 // Determine output range relative to large rectangle
1810 if (lvOutputDimensions.x<cppComplexPlane->iColumns)
1811 {
1812 lrXRange.lMin=0;
1813 lrXRange.lMax=lvOutputDimensions.x;
1814 lXOffset=-lLargeX/2;
1815 }
1816 else
1817 {
1818 lrXRange.lMin=fmax(0, (lLargeX-lvOutputDimensions.x)/2);
1819 lrXRange.lMax=fmin(lvOutputDimensions.x, lrXRange.lMin+lvOutputDimensions.x);
1820 lXOffset=(lLargeX-cppComplexPlane->iColumns)/2;
1821 }
1822 if (lvOutputDimensions.y<cppComplexPlane->iRows)
1823 {
1824 lrYRange.lMin=0;
1825 lrYRange.lMax=lvOutputDimensions.y;
1826 lYOffset=-lLargeY/2;
1827 }
1828 else
1829 {
1830 lrYRange.lMin=fmax(0, (lLargeY-lvOutputDimensions.y)/2);
1831 lrYRange.lMax=fmin(lvOutputDimensions.y, lrYRange.lMin+lvOutputDimensions.y);
1832 lYOffset=(lLargeY-cppComplexPlane->iRows)/2;
1833 }
1834
1835 // Fill output plane
1836 for (i=lrYRange.lMin, y=0; i<lrYRange.lMax; ++i)
1837 {
1838 fpOutputAmplitudeRow=fppOutputAmplitudeBuffer[y];
1839 fpOutputPhaseRow=fppOutputPhaseBuffer[y];
1840
1841 fpInputVector[1] = (float)(i-lHalfRows); // Rotation around centre of image
1842 for (j=lrXRange.lMin, x=0; j<lrXRange.lMax; ++j)
1843 {
1844 fpInputVector[0] = (float)(j-lHalfCols); // Rotation around centre of image
1845
1846 // Determine transformation from point on output image to corresponding point on input image
1847 MatVectMult(mTransMatrix, &(fpInputVector[0]), &(fpOutputVector[0]));
1848
1849 // Correct for offset from UL corner to centre
1850 fpOutputVector[0] += lHalfCols-lXOffset;
1851 fpOutputVector[1] += lHalfRows-lYOffset;
1852
1853 // Find corresponding value if input image point within input image
1854 if (bInterpOrder==NEAREST_NEIGHBOR_INTERPOLATION || fpOutputVector[0]<=iMinX || fpOutputVector[0]>=iMaxX ||
1855 fpOutputVector[1]<=iMinY || fpOutputVector[1]>=iMaxY)
1856 {
1857 FixedPtVector iNearestNeighbor={Round(fpOutputVector[0]), Round(fpOutputVector[1])};
1858
1859 if (iNearestNeighbor.x>=0 && iNearestNeighbor.x<cppComplexPlane->iColumns &&
1860 iNearestNeighbor.y>=0 && iNearestNeighbor.y<cppComplexPlane->iRows)
1861 {
1862 fpOutputAmplitudeRow[j] = (fppAmplitudeBuffer[iNearestNeighbor.y][iNearestNeighbor.x]);
1863 fpOutputPhaseRow[j] = (fppPhaseBuffer[iNearestNeighbor.y][iNearestNeighbor.x]);
1864 }
1865 else
1866 {
1867 if (iNearestNeighbor.x<0) ++(iNearestNeighbor.x);
1868 else if (iNearestNeighbor.x>=cppComplexPlane->iColumns) --(iNearestNeighbor.x);
1869 if (iNearestNeighbor.y<0) ++(iNearestNeighbor.y);
1870 else if (iNearestNeighbor.y>=cppComplexPlane->iRows) --(iNearestNeighbor.y);
1871
1872 if (iNearestNeighbor.x>=0 && iNearestNeighbor.x<cppComplexPlane->iColumns &&
1873 iNearestNeighbor.y>=0 && iNearestNeighbor.y<cppComplexPlane->iRows)
1874 {
1875 fpOutputAmplitudeRow[j] = fppAmplitudeBuffer[iNearestNeighbor.y][iNearestNeighbor.x];
1876 fpOutputPhaseRow[j] = fppPhaseBuffer[iNearestNeighbor.y][iNearestNeighbor.x];
1877 }
1878 else fpOutputAmplitudeRow[j] = fpOutputPhaseRow[j] = 0;
1879 }
1880 }
1881 else
1882 {
1883 fpOutputAmplitudeRow[j] = (
1884 (bInterpOrder==BILINEAR_INTERPOLATION)? (float)Interp_BL(fppAmplitudeBuffer, fpOutputVector[0], fpOutputVector[1]) :
1885 ((bInterpOrder==BIQUADRATIC_INTERPOLATION)? (float)Interp_BQ(fppAmplitudeBuffer, fpOutputVector[0], fpOutputVector[1]) :
1886 (float)Interp_BC(fppAmplitudeBuffer, fpOutputVector[0], fpOutputVector[1])));
1887 fpOutputPhaseRow[j] = (
1888 (bInterpOrder==BILINEAR_INTERPOLATION)? (float)Interp_BL(fppPhaseBuffer, fpOutputVector[0], fpOutputVector[1]) :
1889 ((bInterpOrder==BIQUADRATIC_INTERPOLATION)? (float)Interp_BQ(fppPhaseBuffer, fpOutputVector[0], fpOutputVector[1]) :
1890 (float)Interp_BC(fppPhaseBuffer, fpOutputVector[0], fpOutputVector[1])));
1891 }
1892 ++x;
1893 }
1894 ++y;
1895 }
1896
1897 // Free input plane
1898 FreeComplexPlane(*cppComplexPlane);
1899
1900 // Make new plane with new dimensions
1901 *cppComplexPlane=MakeComplexPlane(lvOutputDimensions.y, lvOutputDimensions.x, NULL);
1902 if (!(cppComplexPlane->Data))
1903 {
1904 for (i=0; i<lvOriginalDimensions.y; ++i)
1905 {
1906 free(fppAmplitudeBuffer[i]);
1907 free(fppPhaseBuffer[i]);
1908 }
1909 for (i=0; i<lvOutputDimensions.y; ++i)
1910 {
1911 free(fppOutputAmplitudeBuffer[i]);
1912 free(fppOutputPhaseBuffer[i]);
1913 }
1914 free(fppAmplitudeBuffer);
1915 free(fppPhaseBuffer);
1916 free(fppOutputAmplitudeBuffer);
1917 free(fppOutputPhaseBuffer);
1918 return ERROR_MEMORY_ALLOCATION;
1919 }
1920 }
1921
1922 // Get complex values from phase and amplitude planes
1923 for (i=0; i<lvOutputDimensions.y; ++i)
1924 {
1925 fpOutputAmplitudeRow=fppOutputAmplitudeBuffer[i];
1926 fpOutputPhaseRow=fppOutputPhaseBuffer[i];
1927 cpOutputRow=cppComplexPlane->Data[i];
1928
1929 for (j=0; j<lvOutputDimensions.x; ++j)
1930 {
1931 /* DEBUG
1932 cpOutputRow[j].real=1;
1933 cpOutputRow[j].imag=1;
1934 */
1935 cpOutputRow[j].real=(float)(fpOutputAmplitudeRow[j]*cos((double)(fpOutputPhaseRow[j])));
1936 cpOutputRow[j].imag=(float)(fpOutputAmplitudeRow[j]*sin((double)(fpOutputPhaseRow[j])));
1937 }
1938 }
1939
1940 for (i=0; i<lvOriginalDimensions.y; ++i)
1941 {
1942 free(fppAmplitudeBuffer[i]);
1943 free(fppPhaseBuffer[i]);
1944 }
1945 for (i=0; i<lvOutputDimensions.y; ++i)
1946 {
1947 free(fppOutputAmplitudeBuffer[i]);
1948 free(fppOutputPhaseBuffer[i]);
1949 }
1950 free(fppAmplitudeBuffer);
1951 free(fppPhaseBuffer);
1952 free(fppOutputAmplitudeBuffer);
1953 free(fppOutputPhaseBuffer);
1954
1955 return ERROR_NONE;
1956 }
1957
1958
MakeCardinalRadialPhaseSamples(FloatPlane fpReal,FloatPlane fpImaginary,short * spCardinalIndices,IntRange lrFrequencyRange,float *** fpppRadialSamples)1959 ERROR_NUMBER MakeCardinalRadialPhaseSamples(FloatPlane fpReal, FloatPlane fpImaginary,
1960 short *spCardinalIndices, IntRange lrFrequencyRange, float ***fpppRadialSamples)
1961 {
1962 ERROR_NUMBER enErrorNumber;
1963 COMPLEX **cpppComplexRadialSample;
1964 int iCenter=fpImaginary.iColumns/2, iRadiusIndex;
1965 int iQuadIndex, iaQuadIndices[4]={iCenter, iCenter, iCenter, iCenter};
1966 size_t stBufferSize=lrFrequencyRange.iMax+1;
1967
1968 if (!(cpppComplexRadialSample=(COMPLEX **)malloc(4*sizeof(COMPLEX *)))) return ERROR_MEMORY_ALLOCATION;
1969 for (iQuadIndex=0; iQuadIndex<4; ++iQuadIndex)
1970 if (!(cpppComplexRadialSample[iQuadIndex]=(COMPLEX *)malloc(stBufferSize*sizeof(COMPLEX))))
1971 {
1972 for (--iQuadIndex; iQuadIndex>=0; --iQuadIndex) free(cpppComplexRadialSample[iQuadIndex]);
1973 free(cpppComplexRadialSample);
1974 return ERROR_MEMORY_ALLOCATION;
1975 }
1976
1977 // Make complex radial samples
1978 for (iRadiusIndex=0; iRadiusIndex<lrFrequencyRange.iMax; ++iRadiusIndex)
1979 {
1980 // 0 degrees (positive along x-axis)
1981 cpppComplexRadialSample[0][iRadiusIndex].real=fpReal.Data[iCenter][iaQuadIndices[0]];
1982 cpppComplexRadialSample[0][iRadiusIndex].imag=fpImaginary.Data[iCenter][iaQuadIndices[0]++];
1983
1984 // 90 degrees (positive along y-axis)
1985 cpppComplexRadialSample[1][iRadiusIndex].real=fpReal.Data[iaQuadIndices[1]][iCenter];
1986 cpppComplexRadialSample[1][iRadiusIndex].imag=fpImaginary.Data[iaQuadIndices[1]++][iCenter];
1987
1988 // 180 degrees (negative along x-axis)
1989 cpppComplexRadialSample[2][iRadiusIndex].real=fpReal.Data[iCenter][iaQuadIndices[2]];
1990 cpppComplexRadialSample[2][iRadiusIndex].imag=fpImaginary.Data[iCenter][iaQuadIndices[2]--];
1991
1992 // 270 degrees (negative along y-axis)
1993 cpppComplexRadialSample[3][iRadiusIndex].real=fpReal.Data[iaQuadIndices[3]][iCenter];
1994 cpppComplexRadialSample[3][iRadiusIndex].imag=fpImaginary.Data[iaQuadIndices[3]--][iCenter];
1995 }
1996
1997 // Make radial phase vector samples from complex radial samples
1998 for (iQuadIndex=0; iQuadIndex<4; ++iQuadIndex)
1999 {
2000 if ((enErrorNumber=MakeRadialPhaseVectorSample(cpppComplexRadialSample[iQuadIndex], lrFrequencyRange.iMax+1,
2001 &((*fpppRadialSamples)[spCardinalIndices[iQuadIndex]])))!=ERROR_NONE)
2002 {
2003 for (iQuadIndex=0; iQuadIndex<4; ++iQuadIndex) free(cpppComplexRadialSample[iQuadIndex]);
2004 free(cpppComplexRadialSample);
2005 return enErrorNumber;
2006 }
2007 }
2008
2009 // Cleanup
2010 for (iQuadIndex=0; iQuadIndex<4; ++iQuadIndex) free(cpppComplexRadialSample[iQuadIndex]);
2011 free(cpppComplexRadialSample);
2012 return ERROR_NONE;
2013 }
2014
2015
ComparePhaseSamples(float * fpTargetSample,float * fpRefSample,IntRange lrFrequencyRange,float * fpComparisonMetric)2016 ERROR_NUMBER ComparePhaseSamples(float *fpTargetSample, float *fpRefSample,
2017 IntRange lrFrequencyRange, float *fpComparisonMetric)
2018 {
2019 ERROR_NUMBER enErrorNumber;
2020 float *fpPhaseShiftVector, fPhaseShiftSlope, fPhaseShiftIntercept, fMSE;
2021 long lVectorLength;
2022
2023 // Make phase shift vector
2024 if ((enErrorNumber=MakePhaseShiftVector(fpTargetSample, fpRefSample, lrFrequencyRange, &lVectorLength,
2025 &fpPhaseShiftVector))!=ERROR_NONE) return enErrorNumber;
2026
2027 // Get slope of phase shift vector
2028 GetSlope(fpPhaseShiftVector, lrFrequencyRange, &fPhaseShiftSlope, &fPhaseShiftIntercept);
2029
2030 // Determine MSE for phase shift vector
2031 fMSE=GetMSE(fpPhaseShiftVector, fPhaseShiftSlope, fPhaseShiftIntercept, lrFrequencyRange);
2032
2033 // Multiplicative inverse of MSE becomes comparison metric
2034 *fpComparisonMetric=(fMSE==0.0f)? 1.0e7f : 1.0f/fMSE;
2035
2036 // Cleanup
2037 free(fpPhaseShiftVector);
2038
2039 return ERROR_NONE;
2040 }
2041 /* Temporarily commented out 2020-09-28
2042 float ExpectedAngularShift(int iCurrentAxis)
2043 {
2044 switch (iCurrentAxis)
2045 {
2046 case X_AXIS:
2047 if (v3DTestRotation.y==0 && v3DTestRotation.z==0) return v3DTestRotation.x;
2048 if (v3DTestRotation.z==0) return v3DTestRotation.x;
2049 return (float)(25.75 * sin(2.0 * DEGREES2RADIANS * v3DTestRotation.x));
2050 case Y_AXIS:
2051 if (v3DTestRotation.x==0 && v3DTestRotation.z==0) return v3DTestRotation.y;
2052 if (v3DTestRotation.z==0) return -v3DTestRotation.y;
2053 return (float)(-90.13 * sin(DEGREES2RADIANS * v3DTestRotation.y));
2054 case Z_AXIS:
2055 if (v3DTestRotation.x==0 && v3DTestRotation.y==0) return v3DTestRotation.z;
2056 if (v3DTestRotation.z==0) return (float)(-v3DTestRotation.z+(24.4*sin(2.0 * DEGREES2RADIANS * v3DTestRotation.z)));
2057 return (float)(25.11 * sin(2.0 * DEGREES2RADIANS * v3DTestRotation.z));
2058 default: return ERROR_PARAMETERS;
2059 }
2060 }
2061 */
2062
GetAngularShift(FloatVector fvCoords)2063 float GetAngularShift(FloatVector fvCoords)
2064 {
2065 float fBuffer;
2066
2067 fBuffer=fvCoords.y-fvCoords.x;
2068 if (fBuffer>180) fBuffer-=360;
2069 else if (fBuffer<-180) fBuffer+=360;
2070
2071 return fBuffer;
2072 }
2073
DetermineUnclippedOutputDimensions(long lInputColumns,long lInputRows,float fClockwiseRotationInDegrees,FloatVector fvTranslation,float fScaleX,float fScaleY)2074 void DetermineUnclippedOutputDimensions(long lInputColumns, long lInputRows,
2075 float fClockwiseRotationInDegrees, FloatVector fvTranslation, float fScaleX, float fScaleY)
2076 {
2077 double dRotationInRadians;
2078 double dCosTheta, dSinTheta;
2079
2080 // Confine rotation angle into (-180,180] degree range.
2081 if (fClockwiseRotationInDegrees>180) fClockwiseRotationInDegrees-=360;
2082 else if (fClockwiseRotationInDegrees<=-180) fClockwiseRotationInDegrees+=360;
2083 dRotationInRadians=DEGREES2RADIANS*fClockwiseRotationInDegrees;
2084
2085 // Determine new dimensions on the basis of rotation
2086 if (fClockwiseRotationInDegrees<-90) // Rotation in (-180,-90) degrees
2087 {
2088 dRotationInRadians=-MATH_PI_2-dRotationInRadians;
2089 dCosTheta=cos(dRotationInRadians);
2090 dSinTheta=sin(dRotationInRadians);
2091 lvOutputDimensions.x=DRound((dSinTheta*lInputColumns)+(dCosTheta*lInputRows)+0.4999999999999);
2092 lvOutputDimensions.y=DRound((dCosTheta*lInputColumns)+(dSinTheta*lInputRows)+0.4999999999999);
2093 }
2094 else if (fClockwiseRotationInDegrees<0) // Rotation in [-90,0) degrees
2095 {
2096 dCosTheta=cos(-dRotationInRadians);
2097 dSinTheta=sin(-dRotationInRadians);
2098 lvOutputDimensions.x=DRound((dCosTheta*lInputColumns)+(dSinTheta*lInputRows)+0.4999999999999);
2099 lvOutputDimensions.y=DRound((dSinTheta*lInputColumns)+(dCosTheta*lInputRows)+0.4999999999999);
2100 }
2101 else if (fClockwiseRotationInDegrees<=90) // Rotation in [0,90] degrees
2102 {
2103 dCosTheta=cos(dRotationInRadians);
2104 dSinTheta=sin(dRotationInRadians);
2105 lvOutputDimensions.x=DRound((dCosTheta*lInputColumns)+(dSinTheta*lInputRows)+0.4999999999999);
2106 lvOutputDimensions.y=DRound((dSinTheta*lInputColumns)+(dCosTheta*lInputRows)+0.4999999999999);
2107 }
2108 else // Rotation in (90,180] degrees
2109 {
2110 dRotationInRadians-=MATH_PI_2;
2111 dCosTheta=cos(dRotationInRadians);
2112 dSinTheta=sin(dRotationInRadians);
2113 lvOutputDimensions.x=DRound((dSinTheta*lInputColumns)+(dCosTheta*lInputRows)+0.4999999999999);
2114 lvOutputDimensions.y=DRound((dCosTheta*lInputColumns)+(dSinTheta*lInputRows)+0.4999999999999);
2115 }
2116
2117 // Determine new dimensions on the basis of scaling
2118 lvOutputDimensions.x=Round(fScaleX*lvOutputDimensions.x);
2119 lvOutputDimensions.y=Round(fScaleY*lvOutputDimensions.y);;
2120
2121 // Determine new dimensions on the basis of translation
2122 lvOutputDimensions.x+=DRound(fabs(fvTranslation.x));
2123 lvOutputDimensions.y+=DRound(fabs(fvTranslation.y));
2124 }
2125
MatVectMult(Matrix mMatrix,float * fpInVector,float * fpOutVector)2126 void MatVectMult(Matrix mMatrix, float *fpInVector, float *fpOutVector)
2127 {
2128 long i, j;
2129 long lMatrixIndex;
2130
2131 for (i=0; i<mMatrix.lRows; ++i)
2132 {
2133 fpOutVector[i] = 0;
2134 lMatrixIndex = i*mMatrix.lColumns;
2135 for (j=0; j<mMatrix.lColumns; ++j)
2136 fpOutVector[i]+=mMatrix.fpData[lMatrixIndex++]*fpInVector[j];
2137 }
2138 }
2139
MakeRadialPhaseVectorSample(COMPLEX * cpComplexRadialSample,int iLength,float ** fppRadialSample)2140 ERROR_NUMBER MakeRadialPhaseVectorSample(COMPLEX *cpComplexRadialSample, int iLength, float **fppRadialSample)
2141 {
2142 float fFrom, fTo, fPhase, fOffset=0, fPrevious=0, fChange;
2143 int iQuad[2], i;
2144 double dThreshold=20, dChangeThresh=MATH_PI_4*3;
2145
2146 // Memory allocation
2147 if (!((*fppRadialSample)=(float *)malloc(iLength*sizeof(float)))) return ERROR_MEMORY_ALLOCATION;
2148 /* DEBUG */
2149
2150 (*fppRadialSample)[0]=0;
2151 for (i=1; i<iLength; ++i)
2152 {
2153 /* DEBUG */
2154
2155 iQuad[1]=(cpComplexRadialSample[i].real>0.0)?
2156 ((cpComplexRadialSample[i].imag>0.0)? 1 : 4) : ((cpComplexRadialSample[i].imag>0.0)? 2 : 3);
2157
2158 fTo=(cpComplexRadialSample[i].real==0)? ((cpComplexRadialSample[i].real>0)? (float)MATH_PI_2 : -(float)MATH_PI_2) :
2159 (float)atan((double)(cpComplexRadialSample[i].imag/cpComplexRadialSample[i].real));
2160 fFrom=fPrevious-fOffset;
2161
2162 if (i>1) switch (iQuad[1])
2163 {
2164 case 1: if (fFrom< -MATH_PI) fTo -= (float)CIRCLE;
2165 if (iQuad[0]==4&&fFrom>MATH_PI)
2166 {
2167 fOffset+=(float)CIRCLE;
2168 }
2169 break;
2170
2171 case 4: if (fFrom> MATH_PI) fTo += (float)CIRCLE;
2172 if (iQuad[0]==1&&fFrom< -(float)MATH_PI)
2173 {
2174 fOffset-=(float)CIRCLE;
2175 }
2176 break;
2177
2178 case 2: if (fFrom> 0) fTo += (float)MATH_PI;
2179 else fTo -= (float)MATH_PI;
2180 break;
2181
2182 case 3: if (fFrom> 0) fTo += (float)MATH_PI;
2183 else fTo -= (float)MATH_PI;
2184 break;
2185 }
2186
2187 fPhase = fOffset+fTo;
2188
2189 // Correct for phase reversal
2190 fChange=fPhase-fPrevious;
2191 if (fChange>dChangeThresh &&
2192 ((FloatPythag((cpComplexRadialSample[i].imag), (cpComplexRadialSample[i].real))<dThreshold) ||
2193 (FloatPythag((cpComplexRadialSample[i-1].imag), (cpComplexRadialSample[i-1].real))<dThreshold)))
2194 {
2195 fOffset-=(float)MATH_PI;
2196 fPhase-=(float)MATH_PI;
2197 }
2198 else if (fChange<-dChangeThresh &&
2199 ((FloatPythag((cpComplexRadialSample[i].imag), (cpComplexRadialSample[i].real))<dThreshold) ||
2200 (FloatPythag((cpComplexRadialSample[i-1].imag), (cpComplexRadialSample[i-1].real))<dThreshold)))
2201 {
2202 fOffset+=(float)MATH_PI;
2203 fPhase+=(float)MATH_PI;
2204 }
2205
2206 iQuad[0]=iQuad[1];
2207 fPrevious=fPhase;
2208 // #if UNIX
2209 #if 1
2210 {
2211 char cDummy;
2212
2213 // Debug
2214 if (fabs((double)fPhase)>1400)
2215 {
2216 fprintf(stderr, "fPhase=%f, i=%d, cpComplexRadialSample[i].fReal=%f\n",
2217 fPhase, i, cpComplexRadialSample[i].real);
2218 fprintf(stderr, "cpComplexRadialSample[i].fImaginary=%f, cpComplexRadialSample[i-1].fReal=%f\n",
2219 cpComplexRadialSample[i].imag, cpComplexRadialSample[i-1].real);
2220 fprintf(stderr, "cpComplexRadialSample[i-1].fImaginary=%f\n", cpComplexRadialSample[i-1].imag);
2221 fscanf(stdin, "%c", &cDummy);
2222 fprintf(stderr, "\n");
2223 }
2224 }
2225 #endif
2226 (*fppRadialSample)[i]=fPhase;
2227 /**/
2228 }
2229 /**/
2230 return ERROR_NONE;
2231 }
2232
MakePhaseShiftVector(float * fpTargetSample,float * fpRefSample,IntRange lrFrequencyRange,long * lpVectorLength,float ** fppPhaseShiftVector)2233 ERROR_NUMBER MakePhaseShiftVector(float *fpTargetSample, float *fpRefSample, IntRange lrFrequencyRange,
2234 long *lpVectorLength, float **fppPhaseShiftVector)
2235 {
2236 ERROR_NUMBER enErrorNumber;
2237 int iInIndex, iOutIndex;
2238 float *fpShiftPtr;
2239 float fPrevious, fDiffUp, fDiffDown, fDiffHere, fOffset=0, fDiffThresh=0.1f;
2240
2241 // Allocate memory to phase shift vector
2242 *lpVectorLength=lrFrequencyRange.iMax+1;
2243 if (!(*fppPhaseShiftVector=(float *)malloc(*lpVectorLength*sizeof(float)))) return ERROR_MEMORY_ALLOCATION;
2244
2245 // Fill vector.
2246 fpShiftPtr=&((*fppPhaseShiftVector)[0]);
2247 for (iOutIndex=0, iInIndex=0; iOutIndex<*lpVectorLength; ++iOutIndex)
2248 {
2249 *fpShiftPtr=fpTargetSample[iInIndex]-fpRefSample[iInIndex]+fOffset;
2250
2251 if (iOutIndex)
2252 {
2253 fDiffUp=(*fpShiftPtr+(float)CIRCLE)-fPrevious;
2254 fDiffDown=fPrevious-*fpShiftPtr+(float)CIRCLE;
2255 fDiffHere=(float)fabs(fPrevious-*fpShiftPtr);
2256
2257 do
2258 {
2259 if (fDiffHere>fDiffDown || fDiffHere>fDiffUp)
2260 {
2261 if (fDiffUp<fDiffDown)
2262 {
2263 *fpShiftPtr += (float)CIRCLE;
2264 fOffset += (float)CIRCLE;
2265 }
2266 else
2267 {
2268 *fpShiftPtr -= (float)CIRCLE;
2269 fOffset -= (float)CIRCLE;
2270 }
2271 }
2272
2273 fDiffUp=(*fpShiftPtr+(float)CIRCLE)-fPrevious;
2274 fDiffDown=fPrevious-*fpShiftPtr+(float)CIRCLE;
2275 fDiffHere=(float)fabs(fPrevious-*fpShiftPtr)-fDiffThresh;
2276 } while (fDiffHere>fDiffDown || fDiffHere>fDiffUp);
2277 }
2278 fPrevious=(iOutIndex<3)? *fpShiftPtr : (2.0f*fPrevious + *fpShiftPtr)/3.0f;
2279 ++fpShiftPtr;
2280 ++iInIndex;
2281 }
2282
2283 // Median filter phase shift vector
2284 if ((enErrorNumber=MedianFilterFloatVector(*fppPhaseShiftVector, *lpVectorLength))!=ERROR_NONE)
2285 {
2286 free(*fppPhaseShiftVector);
2287 return enErrorNumber;
2288 }
2289
2290 return ERROR_NONE;
2291 }
2292
2293
GetSlope(float * fVector,IntRange lrFrequencyRange,float * fpSlope,float * fpIntercept)2294 void GetSlope(float *fVector, IntRange lrFrequencyRange, float *fpSlope, float *fpIntercept)
2295 {
2296 float fXY, fXSquared, *fpShiftPtr;
2297 float fXMean=0, fYMean=0, fXMinusXMean, fYMinusYMean;
2298 int i, iCount=lrFrequencyRange.iMax-lrFrequencyRange.iMin;
2299
2300 // Determine mean X and MeanY
2301 fpShiftPtr=fVector+lrFrequencyRange.iMin;
2302 for (i=lrFrequencyRange.iMin; i<lrFrequencyRange.iMax; ++i, ++fpShiftPtr)
2303 {
2304 fXMean += i;
2305 fYMean += *fpShiftPtr;
2306 }
2307 fXMean /= iCount;
2308 fYMean /= iCount;
2309
2310 fpShiftPtr=fVector+lrFrequencyRange.iMin;
2311 fXY=fXSquared=0.0f;
2312 for (i=lrFrequencyRange.iMin; i<lrFrequencyRange.iMax; ++i, ++fpShiftPtr)
2313 {
2314 fXMinusXMean = (float)i - fXMean;
2315 fYMinusYMean = *fpShiftPtr - fYMean;
2316 fXY += fXMinusXMean * fYMinusYMean;
2317 fXSquared += fXMinusXMean * fXMinusXMean;
2318 }
2319
2320 *fpSlope=fXY/fXSquared;
2321 *fpIntercept=fYMean-(*fpSlope * fXMean);
2322 }
2323
GetMSE(float * fVector,float fSlope,float fPhaseShiftIntercept,IntRange lrFrequencyRange)2324 float GetMSE(float *fVector, float fSlope, float fPhaseShiftIntercept, IntRange lrFrequencyRange)
2325 {
2326 float *fpShiftPtr, fXExpected, fXObserved, fMse, fDiff;
2327 float fOffset=fPhaseShiftIntercept;
2328 int i;
2329
2330 fpShiftPtr=fVector+lrFrequencyRange.iMin;
2331 fMse=0.0;
2332
2333 for (i=lrFrequencyRange.iMin; i<lrFrequencyRange.iMax; ++i, ++fpShiftPtr)
2334 {
2335 fXExpected = (fSlope*i)+fOffset;
2336 fXObserved = *fpShiftPtr;
2337 fDiff=fXObserved-fXExpected;
2338 fMse+=fDiff*fDiff;
2339 }
2340
2341 return(fMse/(float)(lrFrequencyRange.iMax-lrFrequencyRange.iMin));
2342 }
2343
MedianFilterFloatVector(float * fpVector,long lLength)2344 ERROR_NUMBER MedianFilterFloatVector(float *fpVector, long lLength)
2345 {
2346 float *fpBuffer, *fpBufferPtr;
2347 size_t stBytes2Copy=lLength*sizeof(float);
2348 int i, iLast=lLength-1;
2349
2350 // Copy vector into buffer
2351 if (!(fpBuffer=(float *)malloc(stBytes2Copy))) return ERROR_MEMORY_ALLOCATION;
2352 memcpy((void *)fpBuffer, (void *)fpVector, stBytes2Copy);
2353
2354 // Write median filtered version of buffer back into vector.
2355 for (i=1, fpBufferPtr=&(fpBuffer[0]); i<iLast; ++i) fpVector[i]=Float3x1Median(fpBufferPtr++);
2356
2357 // Cleanup
2358 free(fpBuffer);
2359 return ERROR_NONE;
2360 }
2361
Float3x1Median(float * fpInput)2362 float Float3x1Median(float *fpInput)
2363 {
2364 return((fpInput[0]>fpInput[1])? ((fpInput[1]>fpInput[2])? fpInput[1] : (fpInput[0]>fpInput[2])? fpInput[2] : fpInput[0]) :
2365 ((fpInput[1]<fpInput[2])? fpInput[1] : (fpInput[0]<fpInput[2])? fpInput[2] : fpInput[0]));
2366 }
2367
2368
2369 /*******************************************/
FreeFloatPlane(FloatPlane Float_Plane)2370 void FreeFloatPlane(FloatPlane Float_Plane)
2371 /*******************************************/
2372 {
2373 int i;
2374
2375 for (i=0; i<Float_Plane.iRows; ++i) free(Float_Plane.Data[i]);
2376 free(Float_Plane.Data);
2377 Float_Plane.Data = NULL;
2378 }
2379
2380 /*********************/
Round(float arg)2381 int Round(float arg) // Rounds floating point value off to the nearest integer
2382 /*********************/
2383 {
2384 arg+=(arg<0)? -0.5f : 0.5f;
2385 return((int)arg);
2386 }
2387
2388 /***********************************************/
ErrorOpeningFile(char * csFileName)2389 ERROR_NUMBER ErrorOpeningFile(char *csFileName)
2390 /***********************************************/
2391 {
2392 fprintf(stderr, "Error opening %s\n", csFileName);
2393 perror("Reason");
2394
2395 return ERROR_OPENING_FILE;
2396 }
2397
2398 /*******************************************/
FloatPythag(float fArg1,float fArg2)2399 float FloatPythag(float fArg1, float fArg2)
2400 /*******************************************/
2401 {
2402 return (float)(sqrt((double)((fArg1*fArg1)+(fArg2*fArg2))));
2403 }
2404
2405 /******************************************/
Pythag(double dArg1,double dArg2)2406 double Pythag(double dArg1, double dArg2)
2407 /******************************************/
2408 {
2409 return(sqrt((dArg1*dArg1)+(dArg2*dArg2)));
2410 }
2411
2412 /****************************************************************************/
Interp_BC(float ** image,float x_coord,float y_coord)2413 double Interp_BC(float **image, float x_coord, float y_coord)
2414 /****************************************************************************/
2415 /*
2416 Bicubic Interpolation
2417 */
2418 {
2419 int i, j, x[4], y[4], X, Y;
2420
2421 double Phi[4], Psi[4];
2422 double x_x0, x_x1, x_x2, x_x3;
2423 double y_y0, y_y1, y_y2, y_y3;
2424
2425 double Pxy=(double)0;
2426
2427 x[0]=(int)(x_coord+0.5)-1;
2428 y[0]=(int)(y_coord+0.5)-1;
2429 for (i=1; i<4; ++i)
2430 {
2431 x[i]=*x + i;
2432 y[i]=*y + i;
2433 }
2434
2435 x_x0=x_coord-x[0]; /* Numerator factors */
2436 x_x1=x_coord-x[1];
2437 x_x2=x_coord-x[2];
2438 x_x3=x_coord-x[3];
2439
2440 y_y0=y_coord-y[0];
2441 y_y1=y_coord-y[1];
2442 y_y2=y_coord-y[2];
2443 y_y3=y_coord-y[3];
2444
2445 Phi[0]=(x_x1*x_x2*x_x3)/(-6);
2446 Phi[1]=(x_x0*x_x2*x_x3)/(2);
2447 Phi[2]=(x_x0*x_x1*x_x3)/(-2);
2448 Phi[3]=(x_x0*x_x1*x_x2)/(6);
2449
2450 Psi[0]=(y_y1*y_y2*y_y3)/(-6);
2451 Psi[1]=(y_y0*y_y2*y_y3)/(2);
2452 Psi[2]=(y_y0*y_y1*y_y3)/(-2);
2453 Psi[3]=(y_y0*y_y1*y_y2)/(6);
2454
2455 X=x[0]; Y=y[0];
2456
2457 for (i=0; i<4; ++i)
2458 {
2459 for (j=0; j<4; ++j) Pxy += *(*(image+Y)+X++)*Psi[i]*Phi[j];
2460
2461 Y++; X= *x;
2462 }
2463
2464 return(Pxy);
2465 }
2466
2467 /***************************************************************************/
Interp_BQ(float ** image,float x_coord,float y_coord)2468 double Interp_BQ(float **image, float x_coord, float y_coord)
2469 /***************************************************************************/
2470 /*
2471 Biquadratic Interpolation
2472 */
2473 {
2474 int i, j, x[3], y[3], X, Y;
2475
2476 double Phi[3], Psi[3];
2477 double x_x0, x_x1, x_x2;
2478 double y_y0, y_y1, y_y2;
2479
2480 double Pxy=(double)0;
2481
2482 x[0]=(int)(x_coord+0.5)-1;
2483 y[0]=(int)(y_coord+0.5)-1;
2484
2485 for (i=1; i<3; ++i)
2486 {
2487 x[i]=*x + i;
2488 y[i]=*y + i;
2489 }
2490
2491 x_x0=x_coord-x[0]; /* Numerator factors */
2492 x_x1=x_coord-x[1];
2493 x_x2=x_coord-x[2];
2494
2495 y_y0=y_coord-y[0];
2496 y_y1=y_coord-y[1];
2497 y_y2=y_coord-y[2];
2498
2499 Phi[0]=(x_x1*x_x2)/(2);
2500 Phi[1]=(x_x0*x_x2)/(-1);
2501 Phi[2]=(x_x0*x_x1)/(2);
2502
2503 Psi[0]=(y_y1*y_y2)/(2);
2504 Psi[1]=(y_y0*y_y2)/(-1);
2505 Psi[2]=(y_y0*y_y1)/(2);
2506
2507 X=x[0]; Y=y[0];
2508
2509 for (i=0; i<3; ++i)
2510 {
2511 for (j=0; j<3; ++j)
2512 Pxy += *(*(image+Y)+X++)*Psi[i]*Phi[j];
2513
2514 Y++; X= *x;
2515 }
2516
2517 return(Pxy);
2518 }
2519
2520 /**************************************************************/
Interp_BL(float ** image,float x_coord,float y_coord)2521 double Interp_BL(float **image, float x_coord, float y_coord)
2522 /**************************************************************/
2523 /*
2524 Bilinear Interpolation
2525 */
2526 {
2527 int i, j, x[2], y[2], X, Y;
2528
2529 double Phi[2], Psi[2];
2530 double x_x0, x_x1;
2531 double y_y0, y_y1;
2532
2533 double Pxy=(double)0;
2534
2535 x[0]=(int)x_coord;
2536 y[0]=(int)y_coord;
2537
2538 for (i=1; i<2; ++i)
2539 {
2540 x[i]=*x + i;
2541 y[i]=*y + i;
2542 }
2543
2544 x_x0=x_coord-x[0]; /* Numerator factors */
2545 x_x1=x_coord-x[1];
2546
2547 y_y0=y_coord-y[0];
2548 y_y1=y_coord-y[1];
2549
2550 Phi[0]= -x_x1;
2551 Phi[1]=x_x0;
2552
2553 Psi[0]= -y_y1;
2554 Psi[1]=y_y0;
2555
2556 X=x[0]; Y=y[0];
2557
2558 for (i=0; i<2; ++i)
2559 {
2560 for (j=0; j<2; ++j) Pxy += *(*(image+Y)+X++)*Psi[i]*Phi[j];
2561
2562 Y++; X= *x;
2563 }
2564
2565 return(Pxy);
2566 }
2567
2568 /***********************/
DRound(double arg)2569 long DRound(double arg) // Rounds double value off to the nearest integer
2570 /***********************/
2571 {
2572 arg+=(arg<0)? -0.5 : 0.5;
2573 return((long)arg);
2574 }
2575
2576 /*****************************************************************************************************/
Complex2Float(ComplexPlane cpInputPlane,FloatPlane * fppRealOut,FloatPlane * fppImagOut)2577 ERROR_NUMBER Complex2Float(ComplexPlane cpInputPlane, FloatPlane *fppRealOut, FloatPlane *fppImagOut)
2578 /*****************************************************************************************************/
2579 {
2580 int iRow, iCol;
2581 float *fpRealRowPtr, *fpImaginaryRowPtr;
2582 COMPLEX *cpInputRowPtr;
2583
2584 // Make output planes
2585 *fppRealOut=MakeFloatPlane(cpInputPlane.iRows, cpInputPlane.iColumns, NULL);
2586 if (!(fppRealOut->Data)) return ERROR_MEMORY_ALLOCATION;
2587 *fppImagOut=MakeFloatPlane(cpInputPlane.iRows, cpInputPlane.iColumns, NULL);
2588 if (!(fppImagOut->Data))
2589 {
2590 FreeFloatPlane(*fppRealOut);
2591 fppRealOut->Data=NULL;
2592 fppImagOut->Data=NULL;
2593 return ERROR_MEMORY_ALLOCATION;
2594 }
2595
2596 // Fill output planes with data from input plane
2597 for (iRow=0; iRow<cpInputPlane.iRows; ++iRow)
2598 {
2599 fpRealRowPtr=fppRealOut->Data[iRow];
2600 fpImaginaryRowPtr=fppImagOut->Data[iRow];
2601 cpInputRowPtr=cpInputPlane.Data[iRow];
2602 for (iCol=0; iCol<cpInputPlane.iColumns; ++iCol)
2603 {
2604 fpRealRowPtr[iCol]=cpInputRowPtr[iCol].real;
2605 fpImaginaryRowPtr[iCol]=cpInputRowPtr[iCol].imag;
2606 }
2607 }
2608
2609 return ERROR_NONE;
2610 }
2611
2612 /******************************************************************/
MakeFloatPlane(long iRows,long iColumns,float ** Data)2613 FloatPlane MakeFloatPlane(long iRows, long iColumns, float **Data)
2614 /******************************************************************/
2615 {
2616 FloatPlane buffer;
2617 size_t stBytesPerRow=iColumns * sizeof(float);
2618
2619 buffer.iRows=iRows;
2620 buffer.iColumns=iColumns;
2621 buffer.Data=Data;
2622
2623 if ((Data == NULL))
2624 {
2625 if (buffer.Data = (float **)malloc(buffer.iRows * sizeof(float *)))
2626 {
2627 for (long i=0; i<buffer.iRows; ++i)
2628 if ((buffer.Data[i]=(float *)calloc(1, stBytesPerRow)) == NULL)
2629 {
2630 for (--i; i>=0; --i) free(buffer.Data[i]);
2631 free(buffer.Data);
2632 buffer.Data = NULL;
2633 return(buffer);
2634 }
2635 }
2636 }
2637
2638 return(buffer);
2639 }
2640
2641 /**************************************/
RootName(char * csFullPathName)2642 char * RootName(char * csFullPathName)
2643 /**************************************/
2644 {
2645 int iPeriodIndex, iBackslashIndex;
2646 char * FileNameNoExt, * buffer;
2647 int iFullLength=(int)strlen(csFullPathName);
2648
2649 if ((FileNameNoExt=(char *)malloc(iFullLength+1))==NULL) return NULL;
2650 if ((buffer=(char *)malloc(iFullLength+1))==NULL)
2651 {
2652 free(FileNameNoExt);
2653 return NULL;
2654 }
2655 sprintf(FileNameNoExt, "%s", csFullPathName);
2656
2657 if ((iPeriodIndex = (int)(strrchr( csFullPathName, '.' ) - csFullPathName)) > 0)
2658 FileNameNoExt[iPeriodIndex] = '\0';
2659
2660 #if UNIX
2661 if ((iBackslashIndex = strrchr( FileNameNoExt, '/' ) - FileNameNoExt) > 0)
2662 #else
2663 if ((iBackslashIndex = (int)(strrchr( FileNameNoExt, '\\' ) - FileNameNoExt)) > 0)
2664 #endif
2665 sprintf(buffer, "%s", &(FileNameNoExt[iBackslashIndex+1]));
2666 else
2667 sprintf(buffer, "%s", FileNameNoExt);
2668
2669 free(FileNameNoExt);
2670
2671 return (buffer);
2672 }
2673
2674
2675 /*******************************************************************/
GetDirectory(char * csInputFileName,char * csDirectory)2676 ERROR_NUMBER GetDirectory(char *csInputFileName, char *csDirectory)
2677 /*******************************************************************/
2678 {
2679 char *csCurrentDirectory;
2680 long iBackslashIndex;
2681 #if UNIX
2682 csCurrentDirectory;
2683 iBackslashIndex = (long)(strrchr( csInputFileName, '/' ));
2684 #else
2685 char csCurrentDirectory[150];
2686 longlong iBackslashIndex = (longlong)(strrchr( csInputFileName, '\\' ));
2687 #endif
2688
2689 if (iBackslashIndex == 0)
2690 {
2691 #if UNIX
2692 if ((csCurrentDirectory = getcwd(NULL, 150)) == NULL)
2693 {
2694 perror("pwd");
2695 return ERROR_GETTING_DIRECTORY;
2696 }
2697 #else
2698 GetCurrentDirectory( 150, (LPSTR)csCurrentDirectory );
2699 #endif
2700 sprintf(csDirectory, "%s", csCurrentDirectory);
2701 }
2702 else
2703 {
2704 #if UNIX
2705 iBackslashIndex -= (long)csInputFileName;
2706 #else
2707 iBackslashIndex -= (longlong)csInputFileName;
2708 #endif
2709 sprintf(csDirectory, "%s", csInputFileName);
2710 csDirectory[iBackslashIndex] = '\0';
2711 }
2712
2713 return ERROR_NONE;
2714 }
2715