CurvatureScaleSpace

Airplane.txt

64 235 65 234 66 234 67 234 68 233 69 233 70 233 71 232 72 232 73 232 74 232 75 232 76 231 77 231 78 231 79 231 80 231 81 230 82 230 83 230 84 230 85 230 86 230 87 230 88 230 89 230 90 230 91 229 92 229 93 229 94 229 95 229 96 229 97 229 98 229 99 229 100 229 101 228 102 227 101 226 100 226 99 226 99 225 99 224 100 224 101 224 101 225 102 226 103 226 104 226 105 227 106 227 107 227 108 226 109 226 110 226 111 226 112 226 113 225 114 225 115 225 116 225 117 225 118 224 119 224 120 224 121 224 122 224 123 224 124 223 125 223 126 223 127 222 128 222 129 222 130 222 131 222 132 222 133 221 134 221 135 221 136 221 137 221 138 221 139 221 140 221 141 221 142 221 143 221 144 221 145 221 146 221 147 221 148 221 149 221 150 221 151 221 152 221 153 221 154 221 155 221 156 221 157 221 158 221 159 221 160 221 161 221 162 221 163 221 164 221 165 221 166 221 167 221 168 221 169 221 170 221 171 221 172 221 173 221 174 221 175 221 176 221 177 221 178 221 179 221 180 221 181 221 182 221 183 221 184 221 185 221 186 221 187 221 188 221 189 221 190 220 191 220 192 220 193 219 194 219 195 219 196 219 197 219 198 219 199 219 200 219 201 220 202 220 203 220 204 220 205 220 206 220 207 220 208 220 209 220 210 219 211 219 212 219 213 218 214 218 215 217 216 216 217 215 217 214 217 213 217 212 218 211 218 210 218 209 218 208 219 207 220 206 221 205 222 205 223 205 224 204 225 204 226 204 227 204 228 203 229 203 230 203 231 203 232 203 233 203 234 203 235 203 236 203 237 203 238 203 239 203 240 203 241 203 242 203 243 203 244 202 245 202 246 202 247 201 247 200 247 199 248 198 249 198 250 198 251 198 252 197 253 196 254 196 255 195 256 194 257 194 258 193 259 193 260 192 261 192 262 191 263 191 264 190 265 189 266 188 267 188 268 187 269 186 270 186 271 185 272 184 273 184 274 183 275 183 276 182 277 181 278 181 279 180 280 179 281 179 282 178 283 177 284 177 285 176 286 176 287 175 288 174 289 174 290 173 291 172 292 172 293 171 294 171 295 170 296 169 297 169 298 168 299 167 300 167 301 166 302 166 303 165 304 164 305 164 306 163 307 162 308 162 309 161 310 161 311 160 312 159 313 159 314 159 315 158 314 157 313 157 312 157 311 157 310 157 309 157 308 157 307 157 306 158 305 158 304 158 303 158 302 158 301 158 300 158 299 158 298 158 297 158 296 158 295 158 294 158 293 158 292 158 291 158 290 158 289 158 288 158 287 157 286 157 285 157 284 157 283 156 282 156 281 155 280 154 280 153 280 152 280 151 281 150 282 149 283 149 284 148 285 148 286 148 287 147 288 147 289 147 290 147 291 147 292 147 293 146 294 146 295 146 296 146 297 146 298 146 299 146 300 146 301 146 302 146 303 146 304 145 305 145 306 145 307 145 308 145 309 145 310 145 311 145 312 145 313 145 314 145 315 145 316 145 317 145 318 145 319 145 320 145 321 145 322 144 323 144 324 144 325 144 326 144 327 144 328 144 329 144 330 143 330 142 330 141 331 140 331 139 332 138 333 137 334 137 335 137 336 137 337 137 338 136 339 136 340 136 341 136 342 136 343 137 344 137 345 137 346 137 347 136 348 136 349 135 350 135 351 134 352 134 353 133 354 132 355 132 356 131 357 130 358 130 359 129 360 128 361 128 362 127 363 127 364 126 365 125 366 125 367 124 368 123 369 123 370 122 371 122 372 121 373 120 374 120 375 119 376 119 377 118 378 117 379 117 380 116 381 115 382 115 383 114 384 113 385 113 386 112 387 112 388 111 389 110 390 110 391 109 392 109 393 108 394 107 395 107 396 106 397 105 398 105 399 104 400 103 401 103 402 102 403 102 404 101 405 101 406 100 407 99 408 99 409 98 410 97 411 97 412 96 413 96 414 95 415 94 416 94 417 93 418 92 419 92 420 91 421 91 422 90 423 89 424 89 425 88 426 87 426 86 426 85 426 84 426 83 426 82 426 81 427 80 427 79 428 79 429 79 430 79 431 78 432 78 433 78 434 78 435 78 436 78 437 78 438 78 439 77 440 77 441 77 442 77 443 77 444 77 445 77 446 76 447 76 448 75 449 75 450 74 451 74 452 74 453 74 454 74 455 73 456 73 457 73 458 73 459 73 460 73 461 73 462 73 463 74 464 74 465 73 465 72 465 71 466 71 467 71 467 72 467 73 467 74 468 75 469 75 470 76 471 77 471 78 472 79 472 80 472 81 472 82 472 83 472 84 472 85 472 86 472 87 473 88 473 89 473 90 473 91 473 92 473 93 473 94 473 95 473 96 473 97 473 98 473 99 473 100 473 101 473 102 473 103 473 104 473 105 473 106 473 107 473 108 473 109 473 110 473 111 473 112 474 113 474 114 474 115 474 116 474 117 474 118 474 119 474 120 474 121 474 122 474 123 474 124 474 125 474 126 474 127 474 128 474 129 474 130 474 131 474 132 474 133 474 134 474 135 474 136 474 137 474 138 474 139 474 140 474 141 474 142 474 143 474 144 474 145 474 146 474 147 474 148 474 149 474 150 474 151 474 152 474 153 474 154 474 155 474 156 474 157 474 158 474 159 474 160 474 161 474 162 474 163 474 164 474 165 474 166 474 167 474 168 474 169 474 170 474 171 474 172 474 173 474 174 474 175 474 176 474 177 474 178 474 179 474 180 474 181 474 182 474 183 475 184 475 185 475 186 475 187 475 188 475 189 475 190 475 191 475 192 475 193 475 194 475 195 475 196 475 197 475 198 475 199 476 200 476 201 476 202 477 203 478 203 479 203 480 204 481 204 482 203 483 203 484 202 485 202 486 202 487 202 488 202 489 202 490 203 491 203 492 204 493 204 494 205 495 205 496 205 497 206 498 206 499 206 500 207 501 207 502 207 503 208 504 208 505 208 506 209 507 209 508 209 509 209 510 209 511 209 512 209 513 209 514 210 515 210 516 210 517 210 518 211 519 212 520 212 521 212 522 212 523 212 524 212 525 212 526 212 527 212 528 212 529 212 530 212 531 212 532 212 533 213 534 213 535 213 536 213 537 213 538 213 539 213 540 213 541 214 542 214 543 214 544 215 545 216 546 217 547 218 548 219 548 220 549 221 549 222 550 223 550 224 550 225 550 226 551 227 551 228 551 229 551 230 551 231 551 232 551 233 551 234 551 235 551 236 551 237 550 238 550 239 550 240 550 241 549 242 549 243 549 244 549 245 550 246 551 247 552 247 553 248 554 248 555 249 556 249 557 249 558 249 559 249 560 250 561 250 562 250 563 250 564 251 565 251 566 252 567 252 568 253 569 253 570 252 571 252 572 252 573 253 574 254 575 254 576 254 577 255 578 255 579 255 580 256 581 256 582 257 582 258 582 259 582 260 581 260 580 260 579 261 578 262 577 263 576 263 575 264 574 264 573 264 572 264 571 264 570 265 569 265 568 266 567 266 566 266 565 266 564 266 563 266 562 267 561 267 560 267 559 267 558 267 557 267 556 267 555 267 554 267 553 267 552 267 551 267 550 267 549 267 548 267 547 267 546 267 545 267 544 267 543 267 542 267 541 267 540 267 539 267 538 267 537 266 536 266 535 266 534 266 533 266 532 266 531 266 530 266 529 266 528 266 527 265 526 265 525 265 524 264 523 264 522 263 521 263 520 262 519 262 518 262 517 262 516 261 515 261 514 261 513 261 512 261 511 261 510 261 509 261 508 261 507 260 506 260 505 260 504 259 503 259 502 259 501 259 500 259 499 258 498 258 497 258 496 258 495 257 494 257 493 257 492 257 491 257 490 257 489 257 488 257 487 257 486 257 485 257 484 258 483 258 482 258 481 258 480 258 479 258 478 258 477 258 476 258 475 258 474 258 473 259 472 260 472 261 472 262 471 263 471 264 471 265 471 266 471 267 471 268 471 269 471 270 470 271 470 272 470 273 470 274 470 275 470 276 470 277 469 278 469 279 469 280 469 281 469 282 469 283 469 284 469 285 468 286 468 287 468 288 468 289 468 290 468 291 468 292 468 293 467 294 467 295 467 296 467 297 467 298 467 299 467 300 466 301 466 302 466 303 466 304 466 305 466 306 465 307 465 308 465 309 465 310 465 311 465 312 465 313 465 314 464 315 464 316 464 317 464 318 464 319 464 320 464 321 464 322 463 323 463 324 463 325 463 326 463 327 463 328 463 329 462 330 462 331 462 332 462 333 462 334 462 335 462 336 461 337 461 338 461 339 461 340 461 341 461 342 461 343 461 344 460 345 460 346 460 347 460 348 460 349 460 350 459 351 459 352 459 353 459 354 459 355 459 356 458 357 458 358 458 359 458 360 458 361 457 362 457 363 456 364 456 365 455 366 455 367 454 368 453 369 452 370 452 371 451 372 450 372 449 372 448 371 447 371 446 371 445 371 444 371 443 371 442 371 441 371 440 371 439 370 438 370 437 370 436 369 435 369 434 369 433 368 432 368 431 368 430 368 429 368 428 368 427 367 426 367 425 367 424 367 423 367 422 367 421 367 420 367 419 366 418 365 417 364 416 363 416 362 415 361 414 360 413 359 412 359 411 358 410 358 409 357 408 357 407 356 406 356 405 355 404 354 403 354 402 353 401 353 400 352 399 352 398 351 397 351 396 350 395 350 394 349 393 348 392 348 391 348 390 347 389 347 388 346 387 346 386 345 385 344 384 344 383 343 382 343 381 342 380 342 379 341 378 341 377 341 376 340 375 340 374 339 373 338 372 338 371 337 370 337 369 336 368 335 367 335 366 334 365 334 364 333 363 333 362 332 361 332 360 331 359 330 358 329 357 328 356 328 355 327 354 327 353 327 352 326 351 325 350 325 349 324 348 323 347 323 346 322 345 322 344 321 343 320 342 320 341 319 340 319 339 319 338 318 337 317 336 317 335 316 334 316 333 315 332 315 331 314 330 313 329 313 328 312 327 311 326 311 325 311 324 310 323 310 322 310 321 309 320 309 319 309 318 308 317 308 316 308 315 308 314 308 313 309 312 309 311 309 310 309 309 309 308 309 307 309 306 309 305 309 304 309 303 309 302 309 301 309 300 309 299 309 298 309 297 309 296 309 295 309 294 309 293 309 292 309 291 309 290 309 289 309 288 309 287 309 286 309 285 308 284 308 283 308 282 307 281 306 281 305 281 304 281 303 282 302 283 301 284 300 285 300 286 300 287 299 288 299 289 299 290 299 291 299 292 299 293 299 294 299 295 298 296 298 297 298 298 298 299 298 300 298 301 297 300 296 299 296 298 295 297 295 296 294 295 294 294 293 293 293 292 293 291 292 290 291 289 291 288 290 287 290 286 289 285 289 284 288 283 288 282 287 281 286 280 286 279 285 278 285 277 284 276 283 275 283 274 282 273 282 272 281 271 281 270 280 269 280 268 279 267 278 266 277 265 277 264 277 263 276 262 276 261 276 260 275 259 274 258 274 257 273 256 272 255 271 254 271 253 271 252 270 251 270 250 270 249 270 248 270 247 270 246 270 245 270 244 270 243 270 242 270 241 270 240 270 239 270 238 270 237 270 236 270 235 270 234 270 233 270 232 270 231 270 230 270 229 270 228 269 227 269 226 269 225 269 224 269 223 268 222 268 221 267 220 266 219 265 219 264 218 263 218 262 218 261 217 260 217 259 216 258 215 258 214 258 213 259 212 260 211 260 210 260 209 260 208 260 207 260 206 260 205 260 204 260 203 260 202 260 201 260 200 260 199 260 198 260 197 260 196 260 195 260 194 260 193 260 192 260 191 260 190 260 189 260 188 260 187 260 186 260 185 260 184 260 183 260 182 259 181 259 180 259 179 259 178 259 177 258 176 258 175 258 174 258 173 258 172 258 171 258 170 258 169 258 168 258 167 258 166 257 165 257 164 257 163 256 162 256 161 256 160 256 159 256 158 256 157 256 156 255 155 255 154 255 153 255 152 255 151 256 150 256 149 256 148 256 147 255 146 255 145 255 144 255 143 256 142 256 141 256 140 256 139 256 138 256 137 256 136 256 135 255 134 255 133 255 132 255 131 255 130 255 129 255 128 254 127 254 126 254 125 254 124 254 123 254 122 254 121 254 120 254 119 254 118 254 117 254 116 254 115 254 114 253 113 253 112 252 111 251 110 251 109 251 108 250 107 250 106 250 105 250 104 250 103 249 102 248 101 248 100 249 99 249 98 249 97 249 96 249 95 248 94 248 93 248 92 248 91 248 90 248 89 248 88 248 87 247 86 247 85 247 84 247 83 247 82 246 81 246 80 246 79 246 78 246 77 246 76 245 75 245 74 245 73 245 72 244 71 244 70 244 69 243 68 243 67 242 66 242 65 241 65 240 64 239 64 238 64 237 64 236

CurvatureTools.py

"""
Programmer: Chris Tralie
Purpose: To implement curvature/torsion scale space on sampled curves and to 
         create animations of the algorithm in action

Reference:
[1] Mokhtarian, Farzin, and Alan Mackworth. "Scale-based description and recognition of planar curves and two-dimensional shapes." IEEE transactions on pattern analysis and machine intelligence 1 (1986): 34-43.
"""
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from scipy.ndimage.filters import gaussian_filter1d as gf1d
import scipy.io as sio
import scipy.misc
import matplotlib.animation as animation

def getCurvVectors(X, MaxOrder, sigma, loop = False, m = 'nearest'):
    """
    Get smoothed curvature vectors up to a particular order
    Parameters
    ----------
    X: ndarray(N, d)
        An N x d matrix of points in R^d
    MaxOrder: int
        The maximum order of torsion to compute (e.g. 3 for position, velocity, and curvature, and torsion)
    sigma: float
        The smoothing amount
    loop: boolean
        Whether to treat this trajectory as a topological loop (i.e. add an edge between first and last point)
    
    Returns
    -------
    Curvs: A list of (N, 3) arrays, starting with the smoothed curve, then followed
           by the smoothed velocity, curvature, torsion, etc. up to the MaxOrder
    """
    if loop:
        m = 'wrap'
    XSmooth = gf1d(X, sigma, axis=0, order = 0, mode = m)
    Vel = gf1d(X, sigma, axis=0, order = 1, mode = m)
    VelNorm = np.sqrt(np.sum(Vel**2, 1))
    VelNorm[VelNorm == 0] = 1
    Curvs = [XSmooth, Vel]
    for order in range(2, MaxOrder+1):
        Tors = gf1d(X, sigma, axis=0, order=order, mode = m)
        for j in range(1, order):
            #Project away other components
            NormsDenom = np.sum(Curvs[j]**2, 1)
            NormsDenom[NormsDenom == 0] = 1
            Norms = np.sum(Tors*Curvs[j], 1)/NormsDenom
            Tors = Tors - Curvs[j]*Norms[:, None]
        Tors = Tors/(VelNorm[:, None]**order)
        Curvs.append(Tors)
    return Curvs

def getZeroCrossings(Curvs):
    """
    Get zero crossings estimates from all curvature/torsion
    measurements by using the dot product
    Parameters
    ----------
    Curvs: list
    A list of (N, 3) arrays, starting with the smoothed curve, then followed
           by the smoothed velocity, curvature, torsion, etc. up to the MaxOrder
   
    Returns
    -------
    Crossings: list
        List of crossing arrays for each curvature order.  In each array, there
        is a 1 if a sign crossing occurred, zero if otherwise
    """
    Crossings = []
    for C in Curvs:
        dots = np.sum(C[0:-1, :]*C[1::, :], 1)
        cross = np.arange(len(dots))
        cross = cross[dots < 0]
        Crossings.append(cross)
    return Crossings

def getScaleSpaceImages(X, MaxOrder, sigmas, loop):
    """
    Return the curvature scale space images for a sampled spacecurve
    Parameters
    ----------
    X: ndarray(N, d)
        An N x d matrix of points in R^d
    MaxOrder: int
        The maximum order of torsion to compute (e.g. 3 for position, velocity, and curvature, and torsion)
    sigmas: ndarray(M)
        A list of smoothing amounts at which to estimate curvature/torsion/etc
    loop: boolean
        Whether to treat this trajectory as a topological loop (i.e. add an edge between first and last point)
    
    Returns
    -------
    SSImages: list of ndarray(M, N)
        A list of scale space images for each curvature order
    """
    NSigmas = len(sigmas)
    SSImages = []
    for i in range(MaxOrder):
        SSImages.append(np.zeros((NSigmas, X.shape[0])))
    for s in range(len(sigmas)):
        Curvs = getCurvVectors(X, MaxOrder, sigmas[s], loop)
        Crossings = getZeroCrossings(Curvs[1::])
        for i in range(MaxOrder):
            if len(Crossings[i]) > 0:
                SSImages[i][s, Crossings[i]] = 1.0
    return SSImages

class CSSAnimator(animation.FuncAnimation):
    """
    A class for doing animation of curvature scale space images
    for 2D/3D curves
    """
    def __init__(self, fig, X, sigmas, filename, ImageRes = 400, loop = False):
        self.fig = fig
        self.X = X
        self.sigmas = sigmas
        self.ImageRes = ImageRes
        self.SSImage = np.zeros((len(sigmas), X.shape[0]))
        self.loop = loop

        ax1 = fig.add_subplot(221)
        ax2 = fig.add_subplot(222)
        plt.subplot2grid((2, 2), (1, 0), colspan=2)
        ax3 = plt.gca()
        [self.ax1, self.ax2, self.ax3] = [ax1, ax2, ax3]

        #Original curve
        self.origCurve, = ax1.plot(X[:, 0], X[:, 1])
        ax1.set_xticks([])
        ax1.set_yticks([])
        ax1.set_title('Original Curve')

        #Smoothed Curve
        self.smoothCurve = ax2.scatter([0, 0], [0, 0], [1, 1])
        self.smoothCurveInfl = ax2.scatter([], [], 20, 'r')
        ax2.set_xlim([np.min(X[:, 0]), np.max(X[:, 0])])
        ax2.set_ylim([np.min(X[:, 1]), np.max(X[:, 1])])
        plt.title("Smoothed Curve")

        #Scale space image plot
        #I have to hack this so the colormap is scaled properly
        initial = np.zeros((ImageRes, ImageRes))
        initial[0] = 1
        self.ssImagePlot = ax3.imshow(initial, extent = (0, 1, self.sigmas[0], self.sigmas[-1]), interpolation = 'none', aspect = 'auto', cmap=plt.get_cmap('gray'))
        plt.xlabel("t")
        plt.ylabel("Sigma")
        plt.title("Scale Space Image (Zero Crossings)")
        #Setup animation thread
        animation.FuncAnimation.__init__(self, fig, func = self._draw_frame, frames = len(sigmas), interval = 50)

        #Write movie
        FFMpegWriter = animation.writers['ffmpeg']
        metadata = dict(title='Movie Test', artist='Matplotlib',
                        comment='Movie support!')
        writer = FFMpegWriter(fps=15, metadata=metadata, bitrate = 30000)
        self.save("out.mp4", writer = writer)

    def _draw_frame(self, i):
        print i

        Curvs = getCurvVectors(self.X, 2, self.sigmas[i], loop = self.loop)
        Crossings = getZeroCrossings(Curvs)
        XSmooth = Curvs[0]
        Curv = Curvs[2]
        CurvMag = np.sqrt(np.sum(Curv**2, 1)).flatten()
        
        #Draw smoothed curve with inflection points
        self.smoothCurve.set_offsets(XSmooth[:, 0:2])
        self.smoothCurve.set_array(CurvMag)
        self.smoothCurve.__sizes = 20*np.ones(XSmooth.shape[0])
        XInflection = XSmooth[Crossings[2], :]
        self.smoothCurveInfl.set_offsets(XInflection[:, 0:2])
        self.ax2.set_title("Sigma = %.3g"%self.sigmas[i])
        self.ax2.set_xticks([])
        self.ax2.set_yticks([])
        
        self.SSImage[-i-1, Crossings[2]] = 1
        #Do some rudimentary anti-aliasing
        SSImageRes = scipy.misc.imresize(self.SSImage, [self.ImageRes, self.ImageRes])
        SSImageRes[SSImageRes > 0] = 1
        SSImageRes = 1-SSImageRes
        self.ssImagePlot.set_array(SSImageRes)

if __name__ == '__main__':
    X = np.loadtxt('Airplane.txt')
    X = np.reshape(X, (X.size/2, 2))
    Y = np.array(X, dtype=np.float32)
    sigmas = np.linspace(10, 160, 160)

    fig = plt.figure(figsize=(8, 6))
    ani = CSSAnimator(fig, Y, sigmas, "out.mp4", loop = True)