Diff of /APSCpp/APSCpp.c [939879] .. [adbd81]  Maximize  Restore

Switch to unified view

a/APSCpp/APSCpp.c b/APSCpp/APSCpp.c
1
/* APSCpp.cpp  04 Mar 2008 TKS

1
/* APSCpp.cpp  04 Mar 2008 TKS
2
  various routines for enhanced autopano-sift-c

2
  various routines for enhanced autopano-sift-c
3
*/

3
*/
4
/*

4
/*
5
 *  This program is free software; you can redistribute it and/or

5
 *  This program is free software; you can redistribute it and/or
6
 *  modify it under the terms of the GNU General Public

6
 *  modify it under the terms of the GNU General Public
7
 *  License as published by the Free Software Foundation; either

7
 *  License as published by the Free Software Foundation; either
8
 *  version 2 of the License, or (at your option) any later version.

8
 *  version 2 of the License, or (at your option) any later version.
9
 *

9
 *
10
 *  This software is distributed in the hope that it will be useful,

10
 *  This software is distributed in the hope that it will be useful,
11
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of

11
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
12
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU

12
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
13
 *  General Public License for more details.

13
 *  General Public License for more details.
14
 *

14
 *
15
 *  You should have received a copy of the GNU General Public

15
 *  You should have received a copy of the GNU General Public
16
 *  License along with this software; if not, write to the Free Software

16
 *  License along with this software; if not, write to the Free Software
17
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

17
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
18
 *

18
 *
19
 */

19
 */
20
#include "AutoPanoSift.h"

20
#include "AutoPanoSift.h"
21
#include "sphereAlign.h"

21
#include "sphereAlign.h"
22
#include "saInterp.h"

22
#include "saInterp.h"
23
#include <math.h>

23
#include <math.h>
24
#include <tiffio.h>

24
#include <tiffio.h>
25
25
26
// create empty

26
// create empty
27
pDIinfo DIinfo_new0(){

27
pDIinfo DIinfo_new0(){
28
    pDIinfo p = (pDIinfo)malloc(sizeof(DIinfo));

28
    pDIinfo p = (pDIinfo)malloc(sizeof(DIinfo));
29
    if( p ) memset( p, 0, sizeof(DIinfo) );

29
    if( p ) memset( p, 0, sizeof(DIinfo) );
30
    return p;

30
    return p;
31
}

31
}
32
32
33
// delete

33
// delete
34
void DIinfo_delete( pDIinfo p ){

34
void DIinfo_delete( pDIinfo p ){
35
    free( p );

35
    free( p );
36
}

36
}
37
37
38
38
39
/*

39
/*
40
  Get source image info from a panotools project file.

40
  Get source image info from a panotools project file.
41
41
42
  Adapted from "optimizer script parser" in panotools\parser.c

42
  Adapted from "optimizer script parser" in panotools\parser.c
43
  Should work with panotools, ptasmblr and hugin projects.

43
  Should work with panotools, ptasmblr and hugin projects.
44
  Ignores everything except 'i' lines (unless there are none,

44
  Ignores everything except 'i' lines (unless there are none,
45
  when it scans the file again looking for 'o' lines).

45
  when it scans the file again looking for 'o' lines).
46
46
47
  input 

47
  input 
48
     project file name

48
     project file name
49
  output 

49
  output 
50
     adds a DisplayImage struct to DIlist for each "i" or "o" line

50
     adds a DisplayImage struct to DIlist for each "i" or "o" line
51
     with only the following fields set

51
     with only the following fields set
52
       name

52
       name
53
       width

53
       width
54
       height

54
       height
55
       hfov

55
       hfov
56
       format

56
       format
57
       yaw

57
       yaw
58
       pitch

58
       pitch
59
       roll

59
       roll
60
  returns

60
  returns
61
     0 on error else number of image lines found

61
     0 on error else number of image lines found
62
62
63
 Notes

63
 Notes
64
   ignores all optimizable params except hfov, y, p, r; 

64
   ignores all optimizable params except hfov, y, p, r; 
65
   handles back references ("v=n", etc) by copying the value from 

65
   handles back references ("v=n", etc) by copying the value from 
66
   the referenced image.

66
   the referenced image.
67
67
68
   accepts all format codes verbatim (original checks them)

68
   accepts all format codes verbatim (original checks them)
69
*/

69
*/
70
#define LINE_LENGTH 512

70
#define LINE_LENGTH 512
71
#define READ_VAR( fmt, ref )  \

71
#define READ_VAR( fmt, ref )  \
72
    sscanf( ++li, fmt, ref ); \

72
    sscanf( ++li, fmt, ref ); \
73
    if( *li == '"' ) do ++li; while(*li != '"'); \

73
    if( *li == '"' ) do ++li; while(*li != '"'); \
74
    while(*li > ' ') li++;

74
    while(*li > ' ') li++;
75
75
76
int LoadProjectImages( char* script, ArrayList * DIlist )

76
int LoadProjectImages( char* script, ArrayList * DIlist )
77
{

77
{
78
    DIinfo  *im, *ir;  // DIlist items

78
    DIinfo  *im, *ir;  // DIlist items
79
    

79
    
80
    // Variables used by parser

80
    // Variables used by parser
81
    

81
    
82
    char                *li, line[LINE_LENGTH], *ch;

82
    char                *li, line[LINE_LENGTH], *ch;
83
    int                 lineNum = 0;

83
    int                 lineNum = 0;
84
    int                 k;

84
    int                 k;
85
    int                 numIm = 0;

85
    int                 numIm = 0;
86
    int   baseIdx = ArrayList_Count( DIlist );  // 1st image goes here

86
    int   baseIdx = ArrayList_Count( DIlist );  // 1st image goes here
87
87
88
    int keychar;

88
    int keychar;
89
89
90
//    setlocale(LC_ALL, "C");

90
//    setlocale(LC_ALL, "C");
91
// loop over 'i', 'o' linestart

91
// loop over 'i', 'o' linestart
92
    for( keychar = 'i';; ){

92
    for( keychar = 'i';; ){
93
        FILE * fp = fopen( script, "r" );

93
        FILE * fp = fopen( script, "r" );
94
        if( !fp ) return 0;

94
        if( !fp ) return 0;
95
95
96
        // Parse script 

96
        // Parse script 
97
         

97
         
98
        while( (ch = fgets(line, LINE_LENGTH, fp)) )

98
        while( (ch = fgets(line, LINE_LENGTH, fp)) )
99
        {

99
        {
100
            lineNum++;

100
            lineNum++;
101
       

101
       
102
            if( line[0] == keychar 

102
            if( line[0] == keychar 
103
               && isspace( line[1] ) )

103
               && isspace( line[1] ) )
104
            { 

104
            { 
105
            // create new Image descriptor

105
            // create new Image descriptor
106
                im  = DIinfo_new0(); 

106
                im  = DIinfo_new0(); 
107
            // fill it from line      

107
            // fill it from line      
108
                li = &(line[1]);

108
                li = &(line[1]);
109
                while( *li != 0 && *li != '\n')

109
                while( *li != 0 && *li != '\n')
110
                {

110
                {
111
                    switch(*li)

111
                    switch(*li)
112
                    {

112
                    {
113
                        case 'w':   READ_VAR( "%d", &(im->width) )

113
                        case 'w':   READ_VAR( "%d", &(im->width) )
114
                                    break;

114
                                    break;
115
                        case 'h':   READ_VAR( "%d", &(im->height))

115
                        case 'h':   READ_VAR( "%d", &(im->height))
116
                                    break;

116
                                    break;
117
                        case 'v':   if( *(li+1) == '=' ){

117
                        case 'v':   if( *(li+1) == '=' ){
118
                                    // copy referenced value

118
                                    // copy referenced value
119
                                       ++li;

119
                                       ++li;
120
                                       READ_VAR( "%d", &k )

120
                                       READ_VAR( "%d", &k )
121
                                       ir = (DIinfo *)ArrayList_GetItem( DIlist, baseIdx + k );

121
                                       ir = (DIinfo *)ArrayList_GetItem( DIlist, baseIdx + k );
122
                                       im->hfov = ir->hfov;

122
                                       im->hfov = ir->hfov;
123
                                    }else{

123
                                    }else{
124
                                       READ_VAR( "%lf", &(im->hfov));

124
                                       READ_VAR( "%lf", &(im->hfov));
125
                                    }

125
                                    }
126
                                    break;

126
                                    break;
127
                        case 'f':

127
                        case 'f':
128
                                    READ_VAR( "%d", &(im->format) );

128
                                    READ_VAR( "%d", &(im->format) );
129
                                    break;

129
                                    break;
130
                        case 'y':   if( *(li+1) == '=' ){

130
                        case 'y':   if( *(li+1) == '=' ){
131
                                    // copy referenced value

131
                                    // copy referenced value
132
                                       ++li;

132
                                       ++li;
133
                                       READ_VAR( "%d", &k )

133
                                       READ_VAR( "%d", &k )
134
                                       ir = (DIinfo *)ArrayList_GetItem( DIlist, baseIdx + k );

134
                                       ir = (DIinfo *)ArrayList_GetItem( DIlist, baseIdx + k );
135
                                       im->yaw = ir->yaw;

135
                                       im->yaw = ir->yaw;
136
                                    }else{

136
                                    }else{
137
                                        READ_VAR( "%lf", &(im->yaw));

137
                                        READ_VAR( "%lf", &(im->yaw));
138
                                    }

138
                                    }
139
                                    break;

139
                                    break;
140
                        case 'p':   if( *(li+1) == '=' ){

140
                        case 'p':   if( *(li+1) == '=' ){
141
                                    // copy referenced value

141
                                    // copy referenced value
142
                                       ++li;

142
                                       ++li;
143
                                       READ_VAR( "%d", &k )

143
                                       READ_VAR( "%d", &k )
144
                                       ir = (DIinfo *)ArrayList_GetItem( DIlist, baseIdx + k );

144
                                       ir = (DIinfo *)ArrayList_GetItem( DIlist, baseIdx + k );
145
                                       im->pitch = ir->pitch;

145
                                       im->pitch = ir->pitch;
146
                                    }else{

146
                                    }else{
147
                                        READ_VAR("%lf", &(im->pitch));

147
                                        READ_VAR("%lf", &(im->pitch));
148
                                    }

148
                                    }
149
                                    break;

149
                                    break;
150
                        case 'r':   if( *(li+1) == '=' ){

150
                        case 'r':   if( *(li+1) == '=' ){
151
                                    // copy referenced value

151
                                    // copy referenced value
152
                                       ++li;

152
                                       ++li;
153
                                       READ_VAR( "%d", &k )

153
                                       READ_VAR( "%d", &k )
154
                                       ir = (DIinfo *)ArrayList_GetItem( DIlist, baseIdx + k );

154
                                       ir = (DIinfo *)ArrayList_GetItem( DIlist, baseIdx + k );
155
                                       im->roll = ir->roll;

155
                                       im->roll = ir->roll;
156
                                    }else{

156
                                    }else{
157
                                        READ_VAR("%lf", &(im->roll));

157
                                        READ_VAR("%lf", &(im->roll));
158
                                    }

158
                                    }
159
                                    break;

159
                                    break;
160
                        case 'n':           // Set filename

160
                        case 'n':           // Set filename
161
                                {   char *p = im->name;

161
                                {   char *p = im->name;
162
                                    READ_VAR( "%s", p );

162
                                    READ_VAR( "%s", p );
163
                                // strip quotes 

163
                                // strip quotes 
164
                                    if( *p == '"' ){

164
                                    if( *p == '"' ){
165
                                        while( p[1] != '"' ){

165
                                        while( p[1] != '"' ){
166
                                            p[0] = p[1];

166
                                            p[0] = p[1];
167
                                            ++p;

167
                                            ++p;
168
                                        }

168
                                        }
169
                                        *p = 0;

169
                                        *p = 0;
170
                                    }

170
                                    }
171
                                }

171
                                }
172
                                    break;

172
                                    break;
173
                        default:

173
                        default:
174
                            break;

174
                            break;
175
                    }

175
                    }
176
                    while( *li > ' ' ){ // skip next "word"

176
                    while( *li > ' ' ){ // skip next "word"
177
                        if( *li == '"' ) do ++li; while( *li != '"' ); // quoted strings

177
                        if( *li == '"' ) do ++li; while( *li != '"' ); // quoted strings
178
                        ++li;

178
                        ++li;
179
                    }

179
                    }
180
                    if(*li) li++;   // skip one nonnull

180
                    if(*li) li++;   // skip one nonnull
181
181
182
                } // EOL

182
                } // EOL
183
183
184
             // add this image to list

184
             // add this image to list
185
                numIm++;

185
                numIm++;
186
                ArrayList_AddItem( DIlist, im );

186
                ArrayList_AddItem( DIlist, im );
187
187
188
         } // end of line

188
         } // end of line
189
      } // end of file

189
      } // end of file
190
190
191
      fclose( fp );

191
      fclose( fp );
192
192
193
      if( numIm > 0 || keychar == 'o' ) break;

193
      if( numIm > 0 || keychar == 'o' ) break;
194
      keychar = 'o';  // try again...

194
      keychar = 'o';  // try again...
195
  }

195
  }
196
196
197
197
198
  return numIm;

198
  return numIm;
199
199
200
}

200
}
201
201
202
/* Dump an ImagMap to a tiff file

202
/* Dump an ImagMap to a tiff file
203
  return success or failure

203
  return success or failure
204
*/

204
*/
205
int ImageMap_toTIFF( ImageMap * p, char * name )

205
int ImageMap_toTIFF( ImageMap * p, char * name )
206
{

206
{
207
    int i, ok;

207
    int i, ok;
208
208
209
// create tiff file

209
// create tiff file
210
    TIFF * tp = TIFFOpen( name, "w" );

210
    TIFF * tp = TIFFOpen( name, "w" );
211
    if( !tp ) return 0;

211
    if( !tp ) return 0;
212
212
213
// set tiff tags  

213
// set tiff tags  
214
    TIFFSetField( tp, TIFFTAG_IMAGEWIDTH, p->yDim ); 

214
    TIFFSetField( tp, TIFFTAG_IMAGEWIDTH, p->yDim ); 
215
    TIFFSetField( tp, TIFFTAG_IMAGELENGTH, p->xDim ); 

215
    TIFFSetField( tp, TIFFTAG_IMAGELENGTH, p->xDim ); 
216
    TIFFSetField( tp, TIFFTAG_PLANARCONFIG, 1 ); // packed pixels

216
    TIFFSetField( tp, TIFFTAG_PLANARCONFIG, 1 ); // packed pixels
217
    TIFFSetField( tp, TIFFTAG_ROWSPERSTRIP, 1 );  // separate rows

217
    TIFFSetField( tp, TIFFTAG_ROWSPERSTRIP, 1 );  // separate rows
218
    TIFFSetField( tp, TIFFTAG_COMPRESSION, 1 );  // uncompressed

218
    TIFFSetField( tp, TIFFTAG_COMPRESSION, 1 );  // uncompressed
219
    TIFFSetField( tp, TIFFTAG_SAMPLESPERPIXEL, 1 ); // monochrome

219
    TIFFSetField( tp, TIFFTAG_SAMPLESPERPIXEL, 1 ); // monochrome
220
    TIFFSetField( tp, TIFFTAG_PHOTOMETRIC, 1 );     // greyscale

220
    TIFFSetField( tp, TIFFTAG_PHOTOMETRIC, 1 );     // greyscale
221
    TIFFSetField( tp, TIFFTAG_BITSPERSAMPLE, 32 ); 

221
    TIFFSetField( tp, TIFFTAG_BITSPERSAMPLE, 32 ); 
222
    TIFFSetField( tp, TIFFTAG_SAMPLEFORMAT, 3 );    // IEEE float

222
    TIFFSetField( tp, TIFFTAG_SAMPLEFORMAT, 3 );    // IEEE float
223
223
224
// write the data rowwise

224
// write the data rowwise
225
    for( i = 0; i < p->xDim; i++ ){

225
    for( i = 0; i < p->xDim; i++ ){
226
        if((ok = TIFFWriteScanline( tp, p->values[i], i, 0 )) != 1 ) break;

226
        if((ok = TIFFWriteScanline( tp, p->values[i], i, 0 )) != 1 ) break;
227
    }

227
    }
228
228
229
// done

229
// done
230
    TIFFClose( tp );

230
    TIFFClose( tp );
231
    return ok;

231
    return ok;
232
232
233
}

233
}
234
234
235
/*  subroutine equivalent to GenerateKeys.c 

235
/*  subroutine equivalent to GenerateKeys.c 
236
236
237
  for APSC++, new argumnent pdi -> DIinfo struct with the input 

237
  for APSC++, new argumnent pdi -> DIinfo struct with the input 
238
  projection format and angular width into needed to set up a

238
  projection format and angular width into needed to set up a
239
  remapping to stereographic format.  If pdi is 0, behaves like 

239
  remapping to stereographic format.  If pdi is 0, behaves like 
240
  GenerateKeys() with these exceptions.

240
  GenerateKeys() with these exceptions.
241
241
242
  1) it sets the Lowe routine's prefilter sigma to 0.5 if the image stays

242
  1) it sets the Lowe routine's prefilter sigma to 0.5 if the image stays
243
  at original size, or to 0 (no further smoothing) if it is reduced here.

243
  at original size, or to 0 (no further smoothing) if it is reduced here.
244
  

244
  
245
  2) it always reduces image size (when maxdim < largest input dimension) 

245
  2) it always reduces image size (when maxdim < largest input dimension) 
246
  by a smooth interpolation rather than simple decimation.

246
  by a smooth interpolation rather than simple decimation.
247
247
248
  3) it never invokes the Lowe detector's option to expand the image 2x.

248
  3) it never invokes the Lowe detector's option to expand the image 2x.
249
249
250
  returns 0 for failure

250
  returns 0 for failure
251
251
252
  06 Mar 2008 new mode argument.

252
  06 Mar 2008 new mode argument.
253
    0 for legacy (KDTree) mode -- returns list of KeypointN

253
    0 for legacy (KDTree) mode -- returns list of KeypointN
254
    1 for ANN mode -- returns list of Keypoint

254
    1 for ANN mode -- returns list of Keypoint
255
255
256
*/

256
*/
257
257
258
KeypointXMLList * GenerateKeyspp( char * imgname, int maxdim, DIinfo * pdi, int mode )

258
KeypointXMLList * GenerateKeyspp( char * imgname, int maxdim, DIinfo * pdi, int mode )
259
{  

259
{  
260
#ifdef _DEBUG_INTERP

260
#ifdef _DEBUG_INTERP
261
    char tmp[262];

261
    char tmp[262];
262
#endif

262
#endif
263
    int nKeys;

263
    int nKeys;
264
    int pW, pH;

264
    int pW, pH;
265
    double Xcen, Ycen;

265
    double Xcen, Ycen;
266
    saRemap *prm = 0;

266
    saRemap *prm = 0;
267
    ImageMap* picMap = 0;

267
    ImageMap* picMap = 0;
268
    LoweFeatureDetector* lf;

268
    LoweFeatureDetector* lf;
269
    KeypointXMLList* kpp;

269
    KeypointXMLList* kpp;
270
270
271
    double Scale = 1.0;

271
    double Scale = 1.0;
272
272
273
273
274
// reduce Lowe prefilter smoothing

274
// reduce Lowe prefilter smoothing
275
    double savesigma = LoweFeatureDetector_SetPreprocSigma( 0.5 );

275
    double savesigma = LoweFeatureDetector_SetPreprocSigma( 0.5 );
276
276
277
/* read the image */

277
/* read the image */
278
278
279
    DisplayImage* pic = DisplayImage_new(imgname);

279
    DisplayImage* pic = DisplayImage_new(imgname);
280
  // (if we get here the image was read OK, and its name printed)

280
  // (if we get here the image was read OK, and its name printed)
281
    if(!pdi) {

281
    if(!pdi) {
282
        WriteLine("  width %d  height %d", pic->width, pic->height );

282
        WriteLine("  width %d  height %d", pic->width, pic->height );
283
    } else {

283
    } else {
284
        WriteLine("  width %d  height %d  format %d  hfov %g", 

284
        WriteLine("  width %d  height %d  format %d  hfov %g", 
285
            pic->width, pic->height, pdi->format, pdi->hfov );

285
            pic->width, pic->height, pdi->format, pdi->hfov );
286
    // put actual image dimensions in the DIinfo, just in case...

286
    // put actual image dimensions in the DIinfo, just in case...
287
        pdi->width = pic->width; pdi->height = pic->height;

287
        pdi->width = pic->width; pdi->height = pic->height;
288
    }

288
    }
289
289
290
    pW = pic->width;

290
    pW = pic->width;
291
    pH = pic->height;

291
    pH = pic->height;
292
    Xcen = 0.5 * pW;

292
    Xcen = 0.5 * pW;
293
    Ycen = 0.5 * pH;

293
    Ycen = 0.5 * pH;
294
294
295
295
296
// convert to monochrome float and discard input image

296
// convert to monochrome float and discard input image
297
    picMap = DisplayImage_ConvertToImageMap( pic );

297
    picMap = DisplayImage_ConvertToImageMap( pic );
298
    DisplayImage_delete(pic);

298
    DisplayImage_delete(pic);
299
299
300
#ifdef _DEBUG_INTERP

300
#ifdef _DEBUG_INTERP
301
    strcpy(tmp, imgname);

301
    strcpy(tmp, imgname);
302
    strcat(tmp,"-MAP.TIF");

302
    strcat(tmp,"-MAP.TIF");
303
    ImageMap_toTIFF( picMap, tmp );

303
    ImageMap_toTIFF( picMap, tmp );
304
#endif

304
#endif
305
305
306
/* reduce size if required

306
/* reduce size if required
307
  by interpolating in smoothed map

307
  by interpolating in smoothed map
308
  also convert to stereographic if pdi is valid

308
  also convert to stereographic if pdi is valid
309
*/

309
*/
310
  // prescale factor

310
  // prescale factor
311
    if( maxdim > 0) Scale = max( pW, pH ) / (double) maxdim;

311
    if( maxdim > 0) Scale = max( pW, pH ) / (double) maxdim;
312
    if( Scale <= 1 ) Scale = 1;

312
    if( Scale <= 1 ) Scale = 1;
313
313
314
 // interpolate if necessary

314
 // interpolate if necessary
315
    if( pdi || Scale > 1 ){

315
    if( pdi || Scale > 1 ){
316
        ImageMap * pn;

316
        ImageMap * pn;
317
        int x, y;

317
        int x, y;
318
        double X, Y;

318
        double X, Y;
319
        int dwid = pW, dhgt = pH;

319
        int dwid = pW, dhgt = pH;
320
    // reduced dimensions

320
    // reduced dimensions
321
        if( Scale > 1 ){

321
        if( Scale > 1 ){
322
            dwid = (int)(0.5 + pW / Scale );

322
            dwid = (int)(0.5 + pW / Scale );
323
            dhgt = (int)(0.5 + pH / Scale );

323
            dhgt = (int)(0.5 + pH / Scale );
324
            Scale = (double) pW / (double) dwid;    // exact

324
            Scale = (double) pW / (double) dwid;    // exact
325
        WriteLine("  reduce size to %d x %d", dwid, dhgt );

325
        WriteLine("  reduce size to %d x %d", dwid, dhgt );
326
        }

326
        }
327
    // create result map

327
    // create result map
328
        pn = ImageMap_new( dwid, dhgt );

328
        pn = ImageMap_new( dwid, dhgt );
329
    // smooth the source

329
    // smooth the source
330
        //ImageMap_GaussianConvolution( picMap, 0.75 * Scale ); // NO EFFECT.

330
        //ImageMap_GaussianConvolution( picMap, 0.75 * Scale ); // NO EFFECT.
331
        LoweFeatureDetector_SetPreprocSigma( 0 );

331
        LoweFeatureDetector_SetPreprocSigma( 0 );
332
    // set up coordinate mapping

332
    // set up coordinate mapping
333
        if( pdi ){          // ++ mode

333
        if( pdi ){          // ++ mode
334
          CamLens * psp, * pdp;

334
          CamLens * psp, * pdp;
335
        // create the projections

335
        // create the projections
336
          psp = CamLens_new1( pW, pH, pdi->format, pdi->hfov, 16.0 );

336
          psp = CamLens_new1( pW, pH, pdi->format, pdi->hfov, 16.0 );
337
          pdp = CamLens_new1( dwid, dhgt, _stereographic, pdi->hfov, 0 );

337
          pdp = CamLens_new1( dwid, dhgt, _stereographic, pdi->hfov, 0 );
338
        // create the mapping

338
        // create the mapping
339
          prm = saRemap_new( psp, pdp );

339
          prm = saRemap_new( psp, pdp );
340
          saRemap_inv ( prm, 0.5 * dwid, 0.5 * dhgt, &Xcen, &Ycen );

340
          saRemap_inv ( prm, 0.5 * dwid, 0.5 * dhgt, &Xcen, &Ycen );
341
          WriteLine("  use stereographic projection");

341
          WriteLine("  use stereographic projection");
342
        }

342
        }
343
    // interpolate

343
    // interpolate
344
        saInterpSetup((void *)picMap->values,   // void * pdata

344
        saInterpSetup((void *)picMap->values,   // void * pdata
345
                         2,                 // format is column vectors

345
                         2,                 // format is column vectors
346
                         picMap->xDim,      // int wid,

346
                         picMap->xDim,      // int wid,
347
                         picMap->yDim,      // int hgt,

347
                         picMap->yDim,      // int hgt,
348
                         1,                 // samples per pixel

348
                         1,                 // samples per pixel
349
                         picMap->xDim,      // pixels per row

349
                         picMap->xDim,      // pixels per row
350
                         0,                 // char wrap, // required

350
                         0,                 // char wrap, // required
351
                         0,                 // char * pmask,    

351
                         0,                 // char * pmask,    
352
                         0, 0, 0            // mask params

352
                         0, 0, 0            // mask params
353
                      );

353
                      );
354
        for( x = 0; x < pn->xDim; x++ ){

354
        for( x = 0; x < pn->xDim; x++ ){
355
            float * pd = pn->values[x];

355
            float * pd = pn->values[x];
356
            for( y = 0; y < pn->yDim; y++ ){

356
            for( y = 0; y < pn->yDim; y++ ){
357
                X = (double) x; Y = (double) y;     

357
                X = (double) x; Y = (double) y;     
358
                if( prm ) saRemap_inv ( prm, X, Y, &X, &Y );    // ++

358
                if( prm ) saRemap_inv ( prm, X, Y, &X, &Y );    // ++
359
                else { X *= Scale; Y *= Scale; }                // legacy

359
                else { X *= Scale; Y *= Scale; }                // legacy
360
                saInterp_float( &(pd[y]), X, Y );

360
                saInterp_float( &(pd[y]), X, Y );
361
            }

361
            }
362
        }

362
        }
363
    // swap

363
    // swap
364
        ImageMap_delete( picMap );

364
        ImageMap_delete( picMap );
365
        picMap = pn;

365
        picMap = pn;
366
        saRemap_delete( prm ); prm = 0;

366
        saRemap_delete( prm ); prm = 0;
367
#ifdef _DEBUG_INTERP

367
#ifdef _DEBUG_INTERP
368
    strcpy(tmp, imgname);

368
    strcpy(tmp, imgname);
369
    strcat(tmp,"-MAP-PROJ.TIF");

369
    strcat(tmp,"-MAP-PROJ.TIF");
370
    ImageMap_toTIFF( picMap, tmp );

370
    ImageMap_toTIFF( picMap, tmp );
371
    WriteLine("ImageMap dumped to %s", tmp );

371
    WriteLine("ImageMap dumped to %s", tmp );
372
    exit(0);

372
    exit(0);
373
#endif

373
#endif
374
    }

374
    }
375
375
376
376
377
/* find the features 

377
/* find the features 
378
*/

378
*/
379
379
380
    lf = LoweFeatureDetector_new0();

380
    lf = LoweFeatureDetector_new0();
381
    nKeys = LoweFeatureDetector_DetectFeaturesDownscaled (lf, picMap, 0, Scale);

381
    nKeys = LoweFeatureDetector_DetectFeaturesDownscaled (lf, picMap, 0, Scale);
382
    WriteLine("  %d keypoints found\n", nKeys );

382
    WriteLine("  %d keypoints found\n", nKeys );
383
383
384
 /* build the return value

384
 /* build the return value
385
     For KDTree mode, a list of KeypointN; for ANN mode a list of Keypoint, 

385
     For KDTree mode, a list of KeypointN; for ANN mode a list of Keypoint, 
386
     that will be converted to keypointN after matching

386
     that will be converted to keypointN after matching
387
 */

387
 */
388
  {

388
  {
389
    

389
    
390
    ArrayList * globalNaturalKeypoints = 0;

390
    ArrayList * globalNaturalKeypoints = 0;
391
    int i;

391
    int i;
392
392
393
    if( mode == 0 ) globalNaturalKeypoints = ArrayList_new0 ((void *) KeypointN_delete);

393
    if( mode == 0 ) globalNaturalKeypoints = ArrayList_new0 ((void *) KeypointN_delete);
394
394
395
    if( pdi ){ // build remapping function at input scale

395
    if( pdi ){ // build remapping function at input scale
396
      CamLens * psp, * pdp;

396
      CamLens * psp, * pdp;
397
      psp = CamLens_new1( pW, pH, pdi->format, pdi->hfov, 16.0 );

397
      psp = CamLens_new1( pW, pH, pdi->format, pdi->hfov, 16.0 );
398
      pdp = CamLens_new1( pW, pH, _stereographic, pdi->hfov, 0 );

398
      pdp = CamLens_new1( pW, pH, _stereographic, pdi->hfov, 0 );
399
      prm = saRemap_new( psp, pdp );    // dest=>src

399
      prm = saRemap_new( psp, pdp );    // dest=>src
400
    } else prm = 0;

400
    } else prm = 0;
401
401
402
    for(i=0; i < ArrayList_Count(lf->globalKeypoints); i++) {

402
    for(i=0; i < ArrayList_Count(lf->globalKeypoints); i++) {
403
        Keypoint* kp = (Keypoint *) ArrayList_GetItem( lf->globalKeypoints, i );

403
        Keypoint* kp = (Keypoint *) ArrayList_GetItem( lf->globalKeypoints, i );
404
        if( prm ){      // remap coordinates to image projection

404
        if( prm ){      // remap coordinates to image projection
405
            saRemap_inv ( prm,

405
            saRemap_inv ( prm,
406
                          kp->x , kp->y ,   // scale to image size

406
                          kp->x , kp->y ,   // scale to image size
407
                          &kp->x, &kp->y        // remapped coords

407
                          &kp->x, &kp->y        // remapped coords
408
                        );

408
                        );
409
        }

409
        }
410
        if( mode == 0 ) ArrayList_AddItem ( globalNaturalKeypoints, KeypointN_new( kp ));

410
        if( mode == 0 ) ArrayList_AddItem ( globalNaturalKeypoints, KeypointN_new( kp ));
411
    }

411
    }
412
412
413
  // package the list

413
  // package the list
414
    kpp = KeypointXMLList_new0();

414
    kpp = KeypointXMLList_new0();
415
    kpp->imageFile = strdup(imgname);

415
    kpp->imageFile = strdup(imgname);
416
    kpp->xDim = pW;

416
    kpp->xDim = pW;
417
    kpp->yDim = pH;

417
    kpp->yDim = pH;
418
    if( mode == 0 ){

418
    if( mode == 0 ){
419
    // return the KeypointN list

419
    // return the KeypointN list
420
        kpp->array = globalNaturalKeypoints;

420
        kpp->array = globalNaturalKeypoints;
421
    } else {

421
    } else {
422
    // return the Keypoint list

422
    // return the Keypoint list
423
        kpp->array = lf->globalKeypoints;

423
        kpp->array = lf->globalKeypoints;
424
        lf->globalKeypoints = 0;    // don't delete the list!

424
        lf->globalKeypoints = 0;    // don't delete the list!
425
    }

425
    }
426
  }

426
  }
427
  // delete the LoweDetector

427
  // delete the LoweDetector
428
    LoweFeatureDetector_delete(lf);

428
    LoweFeatureDetector_delete(lf);
429
  // restore the default Lowe prefilter

429
  // restore the default Lowe prefilter
430
    LoweFeatureDetector_SetPreprocSigma( savesigma );

430
    LoweFeatureDetector_SetPreprocSigma( savesigma );
431
  // disable the patent warning message for any subsequent images

431
  // disable the patent warning message for any subsequent images
432
    LoweFeatureDetector_SetPrintWarning(false); 

432
    LoweFeatureDetector_SetPrintWarning(false); 
433
433
434
    return kpp;

434
    return kpp;
435
}

435
}
436
436
437
437
438
/* Match keypoints with ANN kd-tree

438
/* Match keypoints with ANN kd-tree
439
  

439
  
440
  argument is a MultiMatch loaded with Keypoints

440
  argument is a MultiMatch loaded with Keypoints
441
  returns with matchlist filled in and the Keypoints

441
  returns with matchlist filled in and the Keypoints
442
  replaced by KeypointN n both keySets and globalKeys

442
  replaced by KeypointN n both keySets and globalKeys
443
  (caution these point to same KeypointN's)

443
  (caution these point to same KeypointN's)
444
444
445
  returns number of matches, or 0 on error

445
  returns number of matches, or 0 on error
446
446
447
*/

447
*/
448
#ifdef __cplusplus

448
#ifdef __cplusplus
449
extern "C" {

449
extern "C" {
450
#endif

450
#endif
451
// wrapper fns in ANNkd_wrap.cpp

451
// wrapper fns in ANNkd_wrap.cpp
452
void ANNkd_new( int sd, double **pts, int n, int d );

452
void ANNkd_new( int sd, double **pts, int n, int d );
453
void ANNkd_search( double * pt, int nnabe, int *indx, double *dist, double nneps );

453
void ANNkd_search( double * pt, int nnabe, int *indx, double *dist, double nneps );
454
void ANNkd_delete( void );

454
void ANNkd_delete( void );
455
#ifdef __cplusplus

455
#ifdef __cplusplus
456
}

456
}
457
#endif

457
#endif
458
458
459
int MultiMatch_ANN ( MultiMatch* self ){

459
int MultiMatch_ANN ( MultiMatch* self ){
460
460
461
// nearest neighbor search parameters

461
// nearest neighbor search parameters
462
#define nnabe 2         // only need 2, self-match is disabled in ANN

462
#define nnabe 2         // only need 2, self-match is disabled in ANN
463
#define nneps 0.0

463
#define nneps 0.0
464
#define mindepth 180    // Nowozin tree uses 130

464
#define mindepth 180    // Nowozin tree uses 130
465
// search results

465
// search results
466
    double dist[nnabe];     // an array of squared distances 

466
    double dist[nnabe];     // an array of squared distances 
467
    int    indx[nnabe];     // an array of point indices

467
    int    indx[nnabe];     // an array of point indices
468
468
469
    int m = 0;  // return value

469
    int m = 0;  // return value
470
    int i, j, k, n, d;

470
    int i, j, k, n, d;
471
    Keypoint * kp;

471
    Keypoint * kp;
472
472
473
  // count keypoints

473
  // count keypoints
474
    n = 0;      // number of keypoints

474
    n = 0;      // number of keypoints
475
    for(i=0; i<ArrayList_Count(self->keySets); i++) {

475
    for(i=0; i<ArrayList_Count(self->keySets); i++) {
476
        KeypointXMLList* list = (KeypointXMLList*) ArrayList_GetItem(self->keySets, i);

476
        KeypointXMLList* list = (KeypointXMLList*) ArrayList_GetItem(self->keySets, i);
477
        n += ArrayList_Count( list->array );

477
        n += ArrayList_Count( list->array );
478
    }

478
    }
479
  /* make array of pointers to coordinates rq'd by ANN

479
  /* make array of pointers to coordinates rq'd by ANN
480
     and array of pointers to corresponding keypointN.s

480
     and array of pointers to corresponding keypointN.s
481
  */

481
  */
482
    double ** pts = (double **)malloc( n * sizeof(double **));

482
    double ** pts = (double **)malloc( n * sizeof(double **));
483
    KeypointN ** kpn = (KeypointN **)malloc( n * sizeof(KeypointN **));

483
    KeypointN ** kpn = (KeypointN **)malloc( n * sizeof(KeypointN **));
484
  // fill pts with pointers to Keypoint feature vectors 

484
  // fill pts with pointers to Keypoint feature vectors 
485
  // and kpn with pointers to equivalent KeypointNs

485
  // and kpn with pointers to equivalent KeypointNs
486
    k = 0;

486
    k = 0;
487
    for(i=0; i<ArrayList_Count(self->keySets); i++) {

487
    for(i=0; i<ArrayList_Count(self->keySets); i++) {
488
        KeypointXMLList* list = (KeypointXMLList*) ArrayList_GetItem(self->keySets, i);

488
        KeypointXMLList* list = (KeypointXMLList*) ArrayList_GetItem(self->keySets, i);
489
        for(j=0; j<ArrayList_Count(list->array); j++) {

489
        for(j=0; j<ArrayList_Count(list->array); j++) {
490
            kp = (Keypoint*) ArrayList_GetItem(list->array, j);

490
            kp = (Keypoint*) ArrayList_GetItem(list->array, j);
491
            pts[k] = kp->featureVector; 

491
            pts[k] = kp->featureVector; 
492
        // Make the KeypointN 

492
        // Make the KeypointN 
493
            kpn[k++] = KeypointN_new( kp );

493
            kpn[k++] = KeypointN_new( kp );
494
        }

494
        }
495
    }

495
    }
496
496
497
    d = kp->featureVectorLength; // # of dimensions

497
    d = kp->featureVectorLength; // # of dimensions
498
498
499
// search cutoff depth 

499
// search cutoff depth 
500
    int sd = (int)( mindepth * max( 1.0, log ((double) n) / log (1000.0) ));

500
    int sd = (int)( mindepth * max( 1.0, log ((double) n) / log (1000.0) ));
501
501
502
  /* build the kd-tree

502
  /* build the kd-tree
503
  */

503
  */
504
504
505
    ANNkd_new( sd, pts, n, d );

505
    ANNkd_new( sd, pts, n, d );
506
    

506
    
507
    if (self->verbose){

507
    if (self->verbose){
508
        WriteLine ("ANN kd-tree: %d keypoints, cutoff depth %d", n, sd);

508
        WriteLine ("ANN kd-tree: %d keypoints, cutoff depth %d", n, sd);
509
    }

509
    }
510
510
511
  // find near neighbors of each key

511
  // find near neighbors of each key
512
    k = 0;

512
    k = 0;
513
    for( i = 0; i < n; i++ ) {

513
    for( i = 0; i < n; i++ ) {
514
    // report progress

514
    // report progress
515
        if (self->verbose) {

515
        if (self->verbose) {
516
            if ((i % 500) == 0)

516
            if ((i % 500) == 0)
517
                Write ("\r%2.2f%%, %d/%d        ",

517
                Write ("\r%2.2f%%, %d/%d        ",
518
                       (100 * ((double) i)) / ((double) n),

518
                       (100 * ((double) i)) / ((double) n),
519
                       i, n);

519
                       i, n);
520
        }

520
        }
521
521
522
    // find neighbors

522
    // find neighbors
523
523
524
        ANNkd_search( pts[i], nnabe, indx, dist, nneps );

524
        ANNkd_search( pts[i], nnabe, indx, dist, nneps );
525
525
526
526
527
    /* add valid neighbor pairs to global matches list

527
    /* add valid neighbor pairs to global matches list
528
        NOTE Matches point to KeypointN's made above

528
        NOTE Matches point to KeypointN's made above
529
        NOTE ANN returns squared distances

529
        NOTE ANN returns squared distances
530
530
531
        Nowozin accepts only the first neighbor as a match, and only if it is 

531
        Nowozin accepts only the first neighbor as a match, and only if it is 
532
        significantly nearer than the 2nd neighbor.  His Match structure includes 

532
        significantly nearer than the 2nd neighbor.  His Match structure includes 
533
        both distances.

533
        both distances.
534
534
535
        The code below disallows matches to points that have already been tested.

535
        The code below disallows matches to points that have already been tested.
536
        That prevents all duplicate matches, at the cost of losing possible control 

536
        That prevents all duplicate matches, at the cost of losing possible control 
537
        points shared by 3 or more images -- but Nowozin's method loses those too.

537
        points shared by 3 or more images -- but Nowozin's method loses those too.
538
    */

538
    */
539
        if( indx[0] > i && dist[0] / dist[1] < 0.36 ){

539
        if( indx[0] > i && dist[0] / dist[1] < 0.36 ){
540
            ArrayList_AddItem( self->globalMatches, 

540
            ArrayList_AddItem( self->globalMatches, 
541
                          Match_new (kpn[i], kpn[indx[0]], 

541
                          Match_new (kpn[i], kpn[indx[0]], 
542
                          sqrt(dist[0]), sqrt(dist[1]))

542
                          sqrt(dist[0]), sqrt(dist[1]))
543
                             );

543
                             );
544
            m++;

544
            m++;
545
        }

545
        }
546
    }

546
    }
547
  

547
  
548
    if (self->verbose) {

548
    if (self->verbose) {
549
        Write ("\r %2.2f%%, %d/%d        ",

549
        Write ("\r %2.2f%%, %d/%d        ",
550
               (100 * ((double) i)) / ((double) n),

550
               (100 * ((double) i)) / ((double) n),
551
               i, n );

551
               i, n );
552
        

552
        
553
        WriteLine ("\nGlobal match search yielded %d matches", m );

553
        WriteLine ("\nGlobal match search yielded %d matches", m );
554
    }

554
    }
555
555
556
  /* Replace Keypoints with KeypointNs in keySets

556
  /* Replace Keypoints with KeypointNs in keySets
557
     NOTE globalKeys list is not needed as only KDTree code uses it

557
     NOTE globalKeys list is not needed as only KDTree code uses it
558
  */

558
  */
559
    k = 0;

559
    k = 0;
560
    for(i=0; i<ArrayList_Count(self->keySets); i++) {

560
    for(i=0; i<ArrayList_Count(self->keySets); i++) {
561
        KeypointXMLList* list = (KeypointXMLList*) ArrayList_GetItem(self->keySets, i);

561
        KeypointXMLList* list = (KeypointXMLList*) ArrayList_GetItem(self->keySets, i);
562
        for(j=0; j<ArrayList_Count(list->array); j++) {

562
        for(j=0; j<ArrayList_Count(list->array); j++) {
563
            kp = (Keypoint*) ArrayList_GetItem(list->array, j);

563
            kp = (Keypoint*) ArrayList_GetItem(list->array, j);
564
            ArrayList_SetItem( list->array, j, (void *)kpn[k] );

564
            ArrayList_SetItem( list->array, j, (void *)kpn[k] );
565
//// Not needed         ArrayList_AddItem( self->globalKeys, (void *)kpn[k] );

565
//// Not needed         ArrayList_AddItem( self->globalKeys, (void *)kpn[k] );
566
            k++;

566
            k++;
567
            Keypoint_delete( kp );

567
            Keypoint_delete( kp );
568
        }

568
        }
569
    }

569
    }
570
570
571
  // tidy

571
  // tidy
572
    free( pts );

572
    free( pts );
573
    free( kpn );

573
    free( kpn );
574
    ANNkd_delete();

574
    ANNkd_delete();
575
575
576
    return m;

576
    return m;
577
}

577
}